Eye size and allometry in deep-sea sergestids
Katie Thomas
February 2021
Data
Morphological data
Morphological measurements (eye diameter and body length) were taken from a mixture of fresh and preserved specimens. In total, 459 sergestid specimens were sampled from 16 species belonging to 12 genera and two subgroups (Sergestes and Sergia). Measurements were recorded to ±0.01mm, but we rounded all to the nearest ±0.1mm as we were not confident in the accuracy of such fine measurements taken with calipers on a rocking ship.
# Full specimen morphology dataset ------
#import data
specimens_raw <- data.frame(read.csv("../Data/specimen_data.csv", header=TRUE, na.strings=c("", "NA", " ", stringsAsFactors = FALSE)))
#tidy data
specimens <- specimens_raw %>%
mutate(Genus = gsub('\\s+', '', Genus)) %>%
mutate(Family = str_to_title(Family)) %>%
mutate(Subgroup = str_to_title(Subgroup)) %>%
mutate(Genus = str_to_title(Genus)) %>%
mutate(genus_species = as.factor(paste(Genus, Species, sep = "_"))) %>%
mutate(pres_bin = recode_factor(Preservation, "ethanol" = "preserved", "paraformaldehyde" = "preserved", "fresh" = "fresh")) %>%
mutate_if(is.numeric,round, 1)
# Species traits ------
#import data
traits_raw <- data.frame(read.csv("../Data/species_data.csv", header=TRUE, na.strings=c("", "NA", " ", stringsAsFactors = FALSE)))
#tidy data
traits <- traits_raw %>%
mutate(Genus = gsub('\\s+', '', Genus)) %>%
mutate(genus_species = as.factor(paste(Genus, Species, sep = "_"))) %>%
mutate(Range = Day_Avg - Night_Avg)
# Show sampling -----
#compile number of shrimp sampled by species
counts <-ddply(specimens, .(specimens$Subgroup, specimens$Genus, specimens$Species), nrow)
#create table column names
names(counts) <- c("Subgroup", "Genus","Species","Sampled")
#generate scrolling table in RMarkdown
kable(counts[ , c("Subgroup", "Genus","Species","Sampled")], caption = "Species and sampling effort for morphological data from family Sergestidae.") %>%
kable_styling(full_width = F) %>%
collapse_rows(columns = 1, valign = "top") %>%
scroll_box(height = "400px") | Subgroup | Genus | Species | Sampled |
|---|---|---|---|
| Sergestes | Allosergestes | pectinatus | 14 |
| Allosergestes | sargassi | 38 | |
| Deosergestes | corniculum | 18 | |
| Deosergestes | henseni | 58 | |
| Eusergestes | arcticus | 5 | |
| Neosergestes | edwardsii | 7 | |
| Parasergestes | armatus | 30 | |
| Parasergestes | vigilax | 17 | |
| Sergestes | atlanticus | 15 | |
| Sergia | Challengerosergia | hansjacobi | 43 |
| Challengerosergia | talismani | 32 | |
| Gardinerosergia | splendens | 42 | |
| Phorcosergia | grandis | 51 | |
| Robustasergia | robusta | 23 | |
| Robustosergia | regalis | 38 | |
| Sergia | tenuiremis | 28 |
Molecular phylogeny
Below is the sergestid molecular phylogeny generated by Charles Golightly. We used ape v.5.4-1 to prune the tree to the 16 species in our dataset and make the tree ultrametric (using the chronopl function with lambda set to 1).
#Import phylogeny
tree_orig <- read.tree(file = "../Data/sertestid_tree_rooted")
#Make copy of tree to manipulate
tree_edit <- tree_orig
#Remove locality tags from phylo tip labels
n <- 3 #which _ to end pattern at
pat <- paste0('^([^_]+(?:_[^_]+){',n-1,'}).*') #pattern
tree_edit$tip.label <- sub(pat, '\\1', tree_edit$tip.label) #remove tags from tips
#Make list of taxa to drop (in tree but not in dataset)
drops <- setdiff(tree_edit$tip.label, as.character(traits$Binomial))
#Drop unwanted tips from phylogeny
tree.pruned <- drop.tip(phy = tree_edit, tip = drops)
#Check and reorder tree
is.rooted(tree.pruned)
is.binary.tree(tree.pruned)
tree <- ladderize(tree.pruned)
#Make row names of the species the phylogeny tip labels
rownames(traits) <- traits$Binomial
#Reorder species data to match phylogeny order
traits <-ReorderData(tree, traits, taxa.names="row names")
#check that phylogeny tips and data match exactly (if they match will return "OK")
name.check(tree, traits)
#make tree ultrametric
tree <- chronopl(tree, 1)#plot tree of sampled species with original tip labels
plot.phylo(tree, #phylogeny to plot
type = "phylogram", #how to visualize phylogeny
show.tip.label = TRUE, #whether to show tip labels/species names
cex = 1.0, #text size
no.margin = TRUE,
use.edge.length = TRUE,
edge.width = 2, #thickness of phylogeny branches
label.offset = 0.1) #how far away from tips to put labels
Static eye-body allometry in Sergestidae
Eye-body allometry in each species
We first examine static post-metamorphosis eye diameter-body length allometry across 16 species of deep-sea sergestids. Allometric relationships were fit using standardised major axis regression via the smatr v.3.4.8 package. For each species, standard model checks, model fits, and an interactive plot of data and fits are shown. Hover over the plots to see individual data values. Points are colored by the source of the sample (fresh collection or preserved).
Note that for 3 of the species, a large/unreasonable outlier was present. We identified these by eye and they have been removed from the data subsets and analyses, but the outlier value is noted under each species header for reference, and a second plot including the outlier as a red diamond is included for transparency (these can go into the supplemental information if we want).
#set colors for fresh and preserved samples
col_pres <- c("ethanol" = "#a6cee3",
"fresh" = "#b2df8a",
"paraformaldehyde" = "#1f78b4")Allosergestes pectinatus
# Subset data -----
apec <- specimens %>% filter(genus_species == "Allosergestes_pectinatus")
# Fit SMA model ------
sma_apec <- sma(formula = Eye_Diameter ~ Body_Length,
data = apec,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_apec, which = "default",type = "o")
plot(sma_apec, which = "residual",type = "o")
plot(sma_apec, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_apec)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = apec, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -2.812398 1.940053
## lower limit -3.491769 1.485663
## upper limit -2.133027 2.533418
##
## H0 : variables uncorrelated
## R-squared : 0.8156794
## P-value : 9.643e-06#save coefficients of fits as object
cc_sma_apec <- data.frame(coef(sma_apec))
# Make plot ------
plot_apec <- ggplot(apec, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Allosergestes pectinatus") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_apec, aes(intercept = cc_sma_apec[1,1], slope = cc_sma_apec[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_apec[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_apec[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_apec)Allosergestes sargassi
Note: this species has a big outlier that was removed for analyses and plots (eye 0.7, body 13.0).
# Subset data -----
#note: outlier removed
asar <- specimens %>%
filter(genus_species == "Allosergestes_sargassi") %>%
filter(!(Eye_Diameter == round(0.67, 1) & Body_Length == round(12.97,1)))
# Fit SMA model ------
sma_asar <- sma(formula = Eye_Diameter ~ Body_Length,
data = asar,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_asar, which = "default",type = "o")
plot(sma_asar, which = "residual",type = "o")
plot(sma_asar, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_asar)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = asar, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -2.728530 1.844597
## lower limit -3.402995 1.419494
## upper limit -2.054065 2.397008
##
## H0 : variables uncorrelated
## R-squared : 0.4037707
## P-value : 2.3852e-05#save coefficients of fits as object
cc_sma_asar <- data.frame(coef(sma_asar))
# Make plot ------
plot_asar <- ggplot(asar, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Allosergestes sargassi") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_asar, aes(intercept = cc_sma_asar[1,1], slope = cc_sma_asar[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_asar[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_asar[1,1], digits = 2), sep = " "), size = 3, col = "black")
#plot with outlier shown
plot_asar_out <- plot_asar +
geom_point(aes(y=round(0.67, 1), x=round(12.97,1)), colour="red", pch = 18) +
ggtitle("Allosergestes sargassi + outlier")
#plot
ggplotly(plot_asar)#plot with outlier
ggplotly(plot_asar_out)Challengerosergia hansjacobi
# Subset data -----
chan <- specimens %>% filter(genus_species == "Challengerosergia_hansjacobi")
# Fit SMA model ------
sma_chan <- sma(formula = Eye_Diameter ~ Body_Length,
data = chan,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_chan, which = "default",type = "o")
plot(sma_chan, which = "residual",type = "o")
plot(sma_chan, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_chan)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = chan, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.528299 1.047036
## lower limit -1.913161 0.824145
## upper limit -1.143436 1.330207
##
## H0 : variables uncorrelated
## R-squared : 0.4129122
## P-value : 3.3885e-06#save coefficients of fits as object
cc_sma_chan <- data.frame(coef(sma_chan))
# Make plot ------
plot_chan <- ggplot(chan, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Challengerosergia hansjacobi") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_chan, aes(intercept = cc_sma_chan[1,1], slope = cc_sma_chan[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_chan[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_chan[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_chan)Challengerosergia talismani
# Subset data -----
ctal <- specimens %>% filter(genus_species == "Challengerosergia_talismani")
# Fit SMA model ------
sma_ctal <- sma(formula = Eye_Diameter ~ Body_Length,
data = ctal,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_ctal, which = "default",type = "o")
plot(sma_ctal, which = "residual",type = "o")
plot(sma_ctal, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_ctal)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = ctal, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.578063 1.0389917
## lower limit -1.840235 0.8845412
## upper limit -1.315891 1.2204109
##
## H0 : variables uncorrelated
## R-squared : 0.8120894
## P-value : 2.0471e-12#save coefficients of fits as object
cc_sma_ctal <- data.frame(coef(sma_ctal))
# Make plot ------
plot_ctal <- ggplot(ctal, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Challengerosergia talismani") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_ctal, aes(intercept = cc_sma_ctal[1,1], slope = cc_sma_ctal[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_ctal[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_ctal[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_ctal)Deosergestes corniculum
# Subset data -----
dcor <- specimens %>% filter(genus_species == "Deosergestes_corniculum")
# Fit SMA model ------
sma_dcor <- sma(formula = Eye_Diameter ~ Body_Length,
data = dcor,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_dcor, which = "default",type = "o")
plot(sma_dcor, which = "residual",type = "o")
plot(sma_dcor, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_dcor)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = dcor, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.432046 0.8653057
## lower limit -1.737583 0.7020666
## upper limit -1.126509 1.0664999
##
## H0 : variables uncorrelated
## R-squared : 0.8421211
## P-value : 8.1768e-08#save coefficients of fits as object
cc_sma_dcor <- data.frame(coef(sma_dcor))
# Make plot ------
plot_dcor <- ggplot(dcor, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Deosergestes corniculum") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_dcor, aes(intercept = cc_sma_dcor[1,1], slope = cc_sma_dcor[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_dcor[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_dcor[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_dcor)Deosergestes henseni
# Subset data -----
dhen <- specimens %>% filter(genus_species == "Deosergestes_henseni")
# Fit SMA model ------
sma_dhen <- sma(formula = Eye_Diameter ~ Body_Length,
data = dhen,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_dhen, which = "default",type = "o")
plot(sma_dhen, which = "residual",type = "o")
plot(sma_dhen, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_dhen)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = dhen, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.306994 0.8238705
## lower limit -1.593941 0.6637531
## upper limit -1.020047 1.0226131
##
## H0 : variables uncorrelated
## R-squared : 0.3380975
## P-value : 1.6982e-06#save coefficients of fits as object
cc_sma_dhen <- data.frame(coef(sma_dhen))
# Make plot ------
plot_dhen <- ggplot(dhen, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Deosergestes henseni") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_dhen, aes(intercept = cc_sma_dhen[1,1], slope = cc_sma_dhen[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_dhen[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_dhen[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_dhen)Eusergestes arcticus
# Subset data -----
euar <- specimens %>% filter(genus_species == "Eusergestes_arcticus")
# Fit SMA model ------
sma_euar <- sma(formula = Eye_Diameter ~ Body_Length,
data = euar,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_euar, which = "default",type = "o")
plot(sma_euar, which = "residual",type = "o")
plot(sma_euar, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_euar)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = euar, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -2.1539322 1.3306114
## lower limit -3.8635100 0.6694485
## upper limit -0.4443543 2.6447543
##
## H0 : variables uncorrelated
## R-squared : 0.8368056
## P-value : 0.029484#save coefficients of fits as object
cc_sma_euar <- data.frame(coef(sma_euar))
# Make plot ------
plot_euar <- ggplot(euar, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Eusergestes arcticus") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_euar, aes(intercept = cc_sma_euar[1,1], slope = cc_sma_euar[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_euar[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_euar[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_euar)Gardinerosergia splendens
Note: this species has a big outlier that was removed for analyses and plots (eye 1.7, body 22.0).
# Subset data -----
gspl <- specimens %>%
filter(genus_species == "Gardinerosergia_splendens") %>%
filter(!(Eye_Diameter == round(1.71, 1) & Body_Length == round(21.98,1)))
# Fit SMA model ------
sma_gspl <- sma(formula = Eye_Diameter ~ Body_Length,
data = gspl,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_gspl, which = "default",type = "o")
plot(sma_gspl, which = "residual",type = "o")
plot(sma_gspl, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_gspl)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = gspl, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.2253950 0.8523703
## lower limit -1.5326713 0.6732695
## upper limit -0.9181186 1.0791149
##
## H0 : variables uncorrelated
## R-squared : 0.4597308
## P-value : 1.1132e-06#save coefficients of fits as object
cc_sma_gspl <- data.frame(coef(sma_gspl))
# Make plot ------
plot_gspl <- ggplot(gspl, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Gardinerosergia splendens") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_gspl, aes(intercept = cc_sma_gspl[1,1], slope = cc_sma_gspl[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_gspl[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_gspl[1,1], digits = 2), sep = " "), size = 3, col = "black")
#plot with outlier shown
plot_gspl_out <- plot_gspl +
geom_point(aes(y=round(1.71,1), x=round(21.98,1)), colour="red", pch = 18) +
ggtitle("Gardinerosergia splendens + outlier")
#plot
ggplotly(plot_gspl)#plot with outlier
ggplotly(plot_gspl_out)Neosergestes edwardsii
Note: this species has a big outlier that was removed for analyses and plots (eye 0.9, body 0.7). Also note that for this species, sample size was insufficient to estimate eye-body allometry (N = 6, p = 0.57 for correlation between eye size and body size).
# Subset data -----
nedw <- specimens %>%
filter(genus_species == "Neosergestes_edwardsii") %>%
filter(!(Eye_Diameter == round(0.85,1) & Body_Length == round(0.72,1)))
# Fit SMA model ------
sma_nedw <- sma(formula = Eye_Diameter ~ Body_Length,
data = nedw,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_nedw, which = "default",type = "o")
plot(sma_nedw, which = "residual",type = "o")
plot(sma_nedw, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_nedw)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = nedw, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate 1.4244363 -1.2712830
## lower limit -0.7533312 -3.8007834
## upper limit 3.6022037 -0.4252177
##
## H0 : variables uncorrelated
## R-squared : 0.08540156
## P-value : 0.57413#save coefficients of fits as object
cc_sma_nedw <- data.frame(coef(sma_nedw))
# Make plot ------
plot_nedw <- ggplot(nedw, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Neosergestes edwardsii") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic"))
#plot with outlier shown
plot_nedw_out <- plot_nedw +
geom_point(aes(x=round(0.85,1), y=round(0.72,1)), colour="red", pch = 18) +
ggtitle("Neosergestes edwardsii + outlier")
#plot
ggplotly(plot_nedw)#plot with outlier
ggplotly(plot_nedw_out)Parasergestes armatus
# Subset data -----
parm <- specimens %>% filter(genus_species == "Parasergestes_armatus")
# Fit SMA model ------
sma_parm <- sma(formula = Eye_Diameter ~ Body_Length,
data = parm,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_parm, which = "default",type = "o")
plot(sma_parm, which = "residual",type = "o")
plot(sma_parm, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_parm)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = parm, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.697577 1.0573117
## lower limit -2.095490 0.8217616
## upper limit -1.299663 1.3603801
##
## H0 : variables uncorrelated
## R-squared : 0.5670657
## P-value : 1.5744e-06#save coefficients of fits as object
cc_sma_parm <- data.frame(coef(sma_parm))
# Make plot ------
plot_parm <- ggplot(parm, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Parasergestes armatus") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_parm, aes(intercept = cc_sma_parm[1,1], slope = cc_sma_parm[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_parm[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_parm[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_parm)Parasergestes vigilax
# Subset data -----
pvig <- specimens %>% filter(genus_species == "Parasergestes_vigilax")
# Fit SMA model ------
sma_pvig <- sma(formula = Eye_Diameter ~ Body_Length,
data = pvig,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_pvig, which = "default",type = "o")
plot(sma_pvig, which = "residual",type = "o")
plot(sma_pvig, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_pvig)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = pvig, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -3.306489 2.375713
## lower limit -4.455793 1.641298
## upper limit -2.157186 3.438750
##
## H0 : variables uncorrelated
## R-squared : 0.5274941
## P-value : 0.00096145#save coefficients of fits as object
cc_sma_pvig <- data.frame(coef(sma_pvig))
# Make plot ------
plot_pvig <- ggplot(pvig, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Parasergestes vigilax") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_pvig, aes(intercept = cc_sma_pvig[1,1], slope = cc_sma_pvig[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_pvig[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_pvig[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_pvig)Phorcosergia grandis
# Subset data -----
pgra <- specimens %>% filter(genus_species == "Phorcosergia_grandis")
# Fit SMA model ------
sma_pgra <- sma(formula = Eye_Diameter ~ Body_Length,
data = pgra,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_pgra, which = "default",type = "o")
plot(sma_pgra, which = "residual",type = "o")
plot(sma_pgra, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_pgra)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = pgra, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.470665 0.9366565
## lower limit -1.632022 0.8480626
## upper limit -1.309308 1.0345055
##
## H0 : variables uncorrelated
## R-squared : 0.8798127
## P-value : < 2.22e-16#save coefficients of fits as object
cc_sma_pgra <- data.frame(coef(sma_pgra))
# Make plot ------
plot_pgra <- ggplot(pgra, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Phorcosergia grandis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_pgra, aes(intercept = cc_sma_pgra[1,1], slope = cc_sma_pgra[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_pgra[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_pgra[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_pgra)Robustasergia robusta
# Subset data -----
rrob <- specimens %>% filter(genus_species == "Robustasergia_robusta")
# Fit SMA model ------
sma_rrob <- sma(formula = Eye_Diameter ~ Body_Length,
data = rrob,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_rrob, which = "default",type = "o")
plot(sma_rrob, which = "residual",type = "o")
plot(sma_rrob, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_rrob)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = rrob, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.2073093 0.8386026
## lower limit -1.5732188 0.6553974
## upper limit -0.8413998 1.0730198
##
## H0 : variables uncorrelated
## R-squared : 0.6989428
## P-value : 6.7833e-07#save coefficients of fits as object
cc_sma_rrob <- data.frame(coef(sma_rrob))
# Make plot ------
plot_rrob <- ggplot(rrob, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Robustasergia robusta") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_rrob, aes(intercept = cc_sma_rrob[1,1], slope = cc_sma_rrob[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_rrob[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_rrob[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_rrob)Robustosergia regalis
# Subset data -----
rreg <- specimens %>% filter(genus_species == "Robustosergia_regalis")
# Fit SMA model ------
sma_rreg <- sma(formula = Eye_Diameter ~ Body_Length,
data = rreg,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_rreg, which = "default",type = "o")
plot(sma_rreg, which = "residual",type = "o")
plot(sma_rreg, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_rreg)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = rreg, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.324440 0.8978007
## lower limit -1.531510 0.7807317
## upper limit -1.117369 1.0324239
##
## H0 : variables uncorrelated
## R-squared : 0.8280327
## P-value : 2.4971e-15#save coefficients of fits as object
cc_sma_rreg <- data.frame(coef(sma_rreg))
# Make plot ------
plot_rreg <- ggplot(rreg, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Robustosergia regalis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_rreg, aes(intercept = cc_sma_rreg[1,1], slope = cc_sma_rreg[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_rreg[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_rreg[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_rreg)Sergestes atlanticus
# Subset data -----
satl <- specimens %>% filter(genus_species == "Sergestes_atlanticus")
# Fit SMA model ------
sma_satl <- sma(formula = Eye_Diameter ~ Body_Length,
data = satl,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_satl, which = "default",type = "o")
plot(sma_satl, which = "residual",type = "o")
plot(sma_satl, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_satl)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = satl, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -2.839127 1.951572
## lower limit -3.727822 1.418514
## upper limit -1.950431 2.684945
##
## H0 : variables uncorrelated
## R-squared : 0.7067609
## P-value : 8.6609e-05#save coefficients of fits as object
cc_sma_satl <- data.frame(coef(sma_satl))
# Make plot ------
plot_satl <- ggplot(satl, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Sergestes atlanticus") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_satl, aes(intercept = cc_sma_satl[1,1], slope = cc_sma_satl[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_satl[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_satl[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_satl)Sergia tenuiremis
# Subset data -----
sten <- specimens %>% filter(genus_species == "Sergia_tenuiremis")
# Fit SMA model ------
sma_sten <- sma(formula = Eye_Diameter ~ Body_Length,
data = sten,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_sten, which = "default",type = "o")
plot(sma_sten, which = "residual",type = "o")
plot(sma_sten, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_sten)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = sten, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.2313110 0.8025570
## lower limit -1.6822103 0.5747311
## upper limit -0.7804117 1.1206941
##
## H0 : variables uncorrelated
## R-squared : 0.2880646
## P-value : 0.0032339#save coefficients of fits as object
cc_sma_sten <- data.frame(coef(sma_sten))
# Make plot ------
plot_sten <- ggplot(sten, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Sergia tenuiremis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_sten, aes(intercept = cc_sma_sten[1,1], slope = cc_sma_sten[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_sten[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_sten[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_sten)Figure for static eye-body scaling across sergestid species
#make figure panels
fig.a <- plot_dcor + theme(legend.position = "none", axis.title.x = element_blank())
fig.b <- plot_dhen + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.c <- plot_apec + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.d <- plot_asar + theme(legend.position = "none", axis.title.x = element_blank())
fig.e <- plot_satl + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.f <- plot_pvig + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.g <- plot_parm + theme(legend.position = "none", axis.title.x = element_blank())
fig.h <- plot_euar + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.i <- plot_gspl + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.j <- plot_rreg + theme(legend.position = "none", axis.title.x = element_blank())
fig.k <- plot_rrob + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.l <- plot_pgra + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.m <- plot_sten + theme(legend.position = "none")
fig.n <- plot_ctal + theme(legend.position = "none", axis.title.y = element_blank())
fig.o <- plot_chan + theme(legend.position = "none", axis.title.y = element_blank())
# arrange panels in figure
plots <- plot_grid(fig.a, fig.b, fig.c, fig.d, fig.e, fig.f, fig.g, fig.h, fig.i, fig.j, fig.k, fig.l, fig.m, fig.n, fig.o, #list of plots to arrange in grid
align = 'vh', #align horizontally and vertically
axis = 'lb',
labels = c("A", "B", "C", "D", "E", "F", "G", "H", "I","J", "K", "L", "M", "N", "O"), #panel labels for figure
hjust = -1, #adjustment for panel labels
nrow = 5) #number of rows in grids
#extract legend
leg <- get_legend(fig.a + theme(legend.position="bottom"))
#view figure
plot_grid(plots, leg, nrow = 2, rel_heights = c(1, .01), align = 'vh')
#export figure
#png("../Figures/dev-allometry.png", width = 10, height = 12, units = "in", res = 500)
pdf("../Figures/static-allometry.pdf", width = 10, height = 12)
plot_grid(plots, leg, nrow = 2, rel_heights = c(1, .05), align = 'vh')
dev.off()Test whether species differ in static allometry
We’ve looked at how eyes scale with body size within species, but are there significant differences in static allometry across species? Here we use two approaches: first, we compare linear mixed models with intercepts varying between species and slopes (+ intercepts) varying between species to see what best describes our data. Then, we look at pairwise comparisons among species to see which are driving these differences.
Linear mixed model to test whether slopes vary across species
Here, we use linear mixed models in the lme4 v.1.1.25 package to test species show differences in static eye-body allometric slopes. We fit 1) a constant-slope, variable intercept model to our eye and body size data with species as a random effect and 2) a variable slope model with the same effects, and then compared model fits by AIC scores. Both models are fit using reduced maximum likelihood (REML), as this approach is better than maximum likelihood (ML) for detecting differences in random effects among models, which is what we are testing here (ML is better for fixed effects).
First we removed the three outliers and the species with insufficient sampling (N. edwardsii) from our full specimen dataset, then ran each model.
#Remove outliers from specimen dataset
specimens_ed <- specimens %>%
filter(!(genus_species == "Allosergestes_sargassi" & Eye_Diameter == 0.7 & Body_Length == 13.0)) %>% #(0.67, 12.97 raw)
filter(!(genus_species == "Gardinerosergia_splendens" & Eye_Diameter == 1.7 & Body_Length == 22.0)) %>% #(1.71, 21.98 raw)
filter(!(genus_species == "Neosergestes_edwardsii" & Eye_Diameter == 0.8 & Body_Length == 0.7)) #(0.85, 0.72 raw)
#Remove species with insufficient sampling
specimens_test <- specimens_ed %>%
filter(genus_species != "Neosergestes_edwardsii") %>%
mutate(species_text = as.factor(paste(Genus, Species, sep = " "))) # Add labels for species textfixed slope model
# variable intercepts model (species can have different intercepts but share mean slope) -------
lme_fixed <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1|genus_species),
data = specimens_test)
# summary
summary(lme_fixed)## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 | genus_species)
## Data: specimens_test
##
## REML criterion at convergence: -1091.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6547 -0.5660 0.0425 0.5699 3.7587
##
## Random effects:
## Groups Name Variance Std.Dev.
## genus_species (Intercept) 0.004257 0.06525
## Residual 0.004549 0.06745
## Number of obs: 450, groups: genus_species, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.24683 0.05121 -24.35
## log10(Body_Length) 0.80642 0.03122 25.83
##
## Correlation of Fixed Effects:
## (Intr)
## lg10(Bdy_L) -0.941# intercepts (variable) and slope (constant)
coef(lme_fixed)## $genus_species
## (Intercept) log10(Body_Length)
## Allosergestes_pectinatus -1.337988 0.8064169
## Allosergestes_sargassi -1.295962 0.8064169
## Challengerosergia_hansjacobi -1.164606 0.8064169
## Challengerosergia_talismani -1.216875 0.8064169
## Deosergestes_corniculum -1.328655 0.8064169
## Deosergestes_henseni -1.278582 0.8064169
## Eusergestes_arcticus -1.246866 0.8064169
## Gardinerosergia_splendens -1.158266 0.8064169
## Parasergestes_armatus -1.325352 0.8064169
## Parasergestes_vigilax -1.299386 0.8064169
## Phorcosergia_grandis -1.246064 0.8064169
## Robustasergia_robusta -1.155267 0.8064169
## Robustosergia_regalis -1.176361 0.8064169
## Sergestes_atlanticus -1.234190 0.8064169
## Sergia_tenuiremis -1.238014 0.8064169
##
## attr(,"class")
## [1] "coef.mer"# AIC
AIC(lme_fixed)## [1] -1083.888#fixed effects
fixef(lme_fixed)## (Intercept) log10(Body_Length)
## -1.2468288 0.8064169#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_fixed, level = 0.95)## 2.5 % 97.5 %
## .sig01 0.04428477 0.09578174
## .sigma 0.06313279 0.07211335
## (Intercept) -1.34880013 -1.14643360
## log10(Body_Length) 0.74536494 0.86929119variable slope model
# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_variable <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species),
data = specimens_test)
# summary
summary(lme_variable)## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |
## genus_species)
## Data: specimens_test
##
## REML criterion at convergence: -1125.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4279 -0.5359 0.0306 0.5910 4.2687
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## genus_species (Intercept) 0.232807 0.48250
## log10(Body_Length) 0.104532 0.32331 -0.99
## Residual 0.003977 0.06306
## Number of obs: 450, groups: genus_species, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.35693 0.13739 -9.876
## log10(Body_Length) 0.89778 0.09181 9.778
##
## Correlation of Fixed Effects:
## (Intr)
## lg10(Bdy_L) -0.993# intercepts (variable) and slope (constant)
coef(lme_variable)## $genus_species
## (Intercept) log10(Body_Length)
## Allosergestes_pectinatus -2.2120964 1.4762026
## Allosergestes_sargassi -1.7513868 1.1366814
## Challengerosergia_hansjacobi -1.0760213 0.7484245
## Challengerosergia_talismani -1.4111004 0.9318641
## Deosergestes_corniculum -1.1953554 0.7273404
## Deosergestes_henseni -0.7843247 0.4963146
## Eusergestes_arcticus -1.2644263 0.8196651
## Gardinerosergia_splendens -0.8716955 0.6169987
## Parasergestes_armatus -1.3258605 0.8075505
## Parasergestes_vigilax -2.0269469 1.3719045
## Phorcosergia_grandis -1.3536721 0.8693323
## Robustasergia_robusta -1.1498050 0.8037740
## Robustosergia_regalis -1.2328226 0.8412590
## Sergestes_atlanticus -1.8945693 1.2766884
## Sergia_tenuiremis -0.8038697 0.5427274
##
## attr(,"class")
## [1] "coef.mer"# AIC
AIC(lme_variable)## [1] -1113.879#fixed effects
fixef(lme_variable)## (Intercept) log10(Body_Length)
## -1.3569302 0.8977818#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_variable, level = 0.95)## 2.5 % 97.5 %
## .sig01 0.27582253 0.74974105
## .sig02 -0.99580190 -0.97348803
## .sig03 0.17701120 0.50987550
## .sigma 0.05903897 0.06767603
## (Intercept) -1.64344107 -1.07983292
## log10(Body_Length) 0.71455335 1.09409513Note that the parameter estimates for slopes and intercepts for species groups are a bit different than those we estimated via OLS. That’s because here, the model is looking across all the data for all the species, and finding a common mean slope/intercept as well as the group effects. Rather than assuming that the species we sampled are all possible species, fitting species as a random effect assumes that the species in our study are some random sampling of the possible species that exist, and that there is some true population mean among species. Here we do see differences across species, but the mixed model pulls the more extreme species toward the group mean (this is all based on variances, both within and between groups/species here). I think this is more appropriate for directly answering the question of whether species are following the same slope for static allometry or not, but it could be argued either way whether these slope and intercept estimates or the ones fitted individually to each species by SMA are more appropriate. If you are trying to predict eye size from body size, I would think that the individual species fit is more appropriate, beause that is a different question. So I think having both methods in the paper is reasonable (each addressing a specific, and different, question), but I am new to mixed models!
compare models
Then we can compare the two models via AIC and log likelihood.
#model AIC comparison
AIC(lme_fixed, lme_variable)## df AIC
## lme_fixed 4 -1083.888
## lme_variable 6 -1113.879#log likelihood comparison
anova(lme_fixed, lme_variable, refit=FALSE)## Data: specimens_test
## Models:
## lme_fixed: log10(Eye_Diameter) ~ log10(Body_Length) + (1 | genus_species)
## lme_variable: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |
## lme_variable: genus_species)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lme_fixed 4 -1083.9 -1067.5 545.94 -1091.9
## lme_variable 6 -1113.9 -1089.2 562.94 -1125.9 33.991 2 4.159e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here we find by comparing AIC scores and log likelihoods that the model with variable slopes and intercepts is better supported than the one with fixed slopes and variable intercepts. We have evidence here that species differ significantly in the slopes of their static allometries rather than following a shared trajectory.
Test for significant effect of preservation as a covariate
Here, we add preservation as a fixed effect into our linear mixed models of eye size vs. body size with species as a random effect. We have already selected the variable slopes model, so here we use a maximum likelihood approach (ML) to fit models that 1) don’t include preservation as a covariate and 2) does include preservation as a covariate and then compare model fits. Note that when we chatted, I mentioned putting in preservation type as a random effect. However, I have read that random effects need to have more than 5 levels, and preservation type has only 3, so by default it needs to be fitted as a fixed effect here.
model without preservation type
# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_nopres <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species),
data = specimens_test,
REML = FALSE) #uses ML instead of REML
# summary
summary(lme_nopres)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |
## genus_species)
## Data: specimens_test
##
## AIC BIC logLik deviance df.resid
## -1123.3 -1098.6 567.6 -1135.3 444
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4474 -0.5397 0.0290 0.5871 4.2627
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## genus_species (Intercept) 0.205896 0.45376
## log10(Body_Length) 0.091875 0.30311 -0.99
## Residual 0.003983 0.06311
## Number of obs: 450, groups: genus_species, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.35193 0.13036 -10.37
## log10(Body_Length) 0.89371 0.08686 10.29
##
## Correlation of Fixed Effects:
## (Intr)
## lg10(Bdy_L) -0.993# intercepts (variable) and slope (constant)
coef(lme_nopres)## $genus_species
## (Intercept) log10(Body_Length)
## Allosergestes_pectinatus -2.1804369 1.4518087
## Allosergestes_sargassi -1.7432805 1.1308540
## Challengerosergia_hansjacobi -1.0817261 0.7520710
## Challengerosergia_talismani -1.4101468 0.9312295
## Deosergestes_corniculum -1.1915659 0.7253709
## Deosergestes_henseni -0.7877675 0.4985130
## Eusergestes_arcticus -1.2643293 0.8197731
## Gardinerosergia_splendens -0.8765113 0.6200664
## Parasergestes_armatus -1.3286959 0.8097119
## Parasergestes_vigilax -1.9939820 1.3460557
## Phorcosergia_grandis -1.3527444 0.8688290
## Robustasergia_robusta -1.1602353 0.8095809
## Robustosergia_regalis -1.2356690 0.8429100
## Sergestes_atlanticus -1.8546819 1.2481093
## Sergia_tenuiremis -0.8172020 0.5507936
##
## attr(,"class")
## [1] "coef.mer"# AIC
AIC(lme_nopres)## [1] -1123.278#fixed effects
fixef(lme_nopres)## (Intercept) log10(Body_Length)
## -1.3519317 0.8937118#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_nopres, level = 0.95)## 2.5 % 97.5 %
## .sig01 0.27582245 0.74974090
## .sig02 -0.99580205 -0.97348803
## .sig03 0.17701116 0.50987543
## .sigma 0.05903897 0.06767603
## (Intercept) -1.64344108 -1.07983293
## log10(Body_Length) 0.71455335 1.09409513model including preservation type
#reorder levels for preservation so that contrasts will be compared to fresh samples
specimens_test$Preservation <- factor(specimens_test$Preservation,
levels = c("fresh", "ethanol", "paraformaldehyde"))
# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_withpres <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 + log10(Body_Length) | genus_species),
data = specimens_test,
REML = FALSE) #uses ML instead of REML
# summary
summary(lme_withpres)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 +
## log10(Body_Length) | genus_species)
## Data: specimens_test
##
## AIC BIC logLik deviance df.resid
## -1152.9 -1120.0 584.5 -1168.9 442
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7037 -0.5122 0.0295 0.5644 4.3676
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## genus_species (Intercept) 0.167876 0.40973
## log10(Body_Length) 0.071632 0.26764 -0.99
## Residual 0.003701 0.06083
## Number of obs: 450, groups: genus_species, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.353719 0.119794 -11.300
## log10(Body_Length) 0.894613 0.078169 11.445
## Preservationethanol -0.019607 0.008151 -2.405
## Preservationparaformaldehyde 0.045520 0.012167 3.741
##
## Correlation of Fixed Effects:
## (Intr) l10(B_ Prsrvtnt
## lg10(Bdy_L) -0.990
## Prsrvtnthnl -0.057 0.022
## Prsrvtnprfr -0.068 0.041 0.455# AIC
AIC(lme_withpres)## [1] -1152.918#fixed effects
fixef(lme_withpres)## (Intercept) log10(Body_Length)
## -1.35371863 0.89461257
## Preservationethanol Preservationparaformaldehyde
## -0.01960702 0.04551984#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_withpres, level = 0.95)## 2.5 % 97.5 %
## .sig01 0.25281798 0.67323346
## .sig02 -0.99341603 -0.96585906
## .sig03 0.15905142 0.44860338
## .sigma 0.05692288 0.06521125
## (Intercept) -1.61782966 -1.10480110
## log10(Body_Length) 0.73401476 1.07261072
## Preservationethanol -0.03573800 -0.00350345
## Preservationparaformaldehyde 0.02141721 0.06962976compare models
Then we can compare the two models via AIC and log likelihood.
#model AIC comparison
AIC(lme_nopres, lme_withpres)## df AIC
## lme_nopres 6 -1123.278
## lme_withpres 8 -1152.918#log likelihood comparison
anova(lme_nopres, lme_withpres)## Data: specimens_test
## Models:
## lme_nopres: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |
## lme_nopres: genus_species)
## lme_withpres: log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 +
## lme_withpres: log10(Body_Length) | genus_species)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lme_nopres 6 -1123.3 -1098.6 567.64 -1135.3
## lme_withpres 8 -1152.9 -1120.0 584.46 -1168.9 33.64 2 4.957e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here we see that the model that includes preservation type as a covariate is significantly better supported. If you take a look at the parameter estimates from the model, you can see that ethanol preservation has a negative effect on intercept (eye size compared to body size), but paraformaldehyde has a positive effect. Note also that these effects are very small, so while they are significantly affecting our models, they seem to have minor effects. I would use this to argue why we don’t include preservation type in all our ecology ancovas; our power is already low due to small sample sizes and including preservation as a covariate would probalby overparamaterize our models.
Linear SMA regression for pairwise comparisons among species slopes
The smatr package also allows you to test whether groups exhibit significant differences in allometric scaling (slopes or intercepts). Below, we again use standardized major axis regression to model the scaling of eye diameter with body length across all specimens measured and test whether we see significant differences in ontogenetic scaling across species.
*Note: As we did not have sufficient sampling to estimate eye-body allometry for Neosergestes edwardsii, I’ve removed it from the specimen dataset before fitting this model.
This test shows that among 15 sergestid species (450 specimens), species significantly differ in their allometric scaling relationships between eye diameter and body length (p < 0.0001). It also yields a table of pairwise comparisons so you can see exactly which species differ significantly in slope and which do not.
# Model for static allometry across specimens with species and species*body size interaction as covariates
sma_species <- sma(formula = Eye_Diameter ~ Body_Length * genus_species,
data = specimens_test,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
multcomp = TRUE, # adds pairwise comparisons between groups
multcompmethod = "adjusted", # adjusts p-value for multiple comparisons
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_species, which = "default",type = "o")
plot(sma_species,which = "residual",type = "o")
plot(sma_species, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_species)## Call: sma(formula = Eye_Diameter ~ Body_Length * genus_species, data = specimens_test,
## log = "xy", method = "SMA", alpha = 0.05, multcomp = TRUE,
## multcompmethod = "adjusted")
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Results of comparing lines among groups.
##
## H0 : slopes are equal.
## Likelihood ratio statistic : 78.37 with 14 degrees of freedom
## P-value : 5.6738e-11
## ------------------------------------------------------------
##
## Results of multiple comparisons among groups.
##
## Test for pair-wise difference in slope :
## genus_species_1 genus_species_2 Pval
## 1 Allosergestes_pectinatus Allosergestes_sargassi 1.00000000
## 2 Allosergestes_pectinatus Challengerosergia_hansjacobi 0.10631965
## 3 Allosergestes_pectinatus Challengerosergia_talismani 0.03771725
## 4 Allosergestes_pectinatus Deosergestes_corniculum 0.00370375
## 5 Allosergestes_pectinatus Deosergestes_henseni 0.00191135
## 6 Allosergestes_pectinatus Eusergestes_arcticus 1.00000000
## 7 Allosergestes_pectinatus Gardinerosergia_splendens 0.00417512
## 8 Allosergestes_pectinatus Parasergestes_armatus 0.14662731
## 9 Allosergestes_pectinatus Parasergestes_vigilax 1.00000000
## 10 Allosergestes_pectinatus Phorcosergia_grandis 0.00484593
## 11 Allosergestes_pectinatus Robustasergia_robusta 0.00386781
## 12 Allosergestes_pectinatus Robustosergia_regalis 0.00330606
## 13 Allosergestes_pectinatus Sergestes_atlanticus 1.00000000
## 14 Allosergestes_pectinatus Sergia_tenuiremis 0.00977338
## 15 Allosergestes_sargassi Challengerosergia_hansjacobi 0.17642048
## 16 Allosergestes_sargassi Challengerosergia_talismani 0.03323729
## 17 Allosergestes_sargassi Deosergestes_corniculum 0.00219121
## 18 Allosergestes_sargassi Deosergestes_henseni 0.00062265
## 19 Allosergestes_sargassi Eusergestes_arcticus 1.00000000
## 20 Allosergestes_sargassi Gardinerosergia_splendens 0.00279183
## 21 Allosergestes_sargassi Parasergestes_armatus 0.25381226
## 22 Allosergestes_sargassi Parasergestes_vigilax 1.00000000
## 23 Allosergestes_sargassi Phorcosergia_grandis 0.00071873
## 24 Allosergestes_sargassi Robustasergia_robusta 0.00340588
## 25 Allosergestes_sargassi Robustosergia_regalis 0.00055799
## 26 Allosergestes_sargassi Sergestes_atlanticus 1.00000000
## 27 Allosergestes_sargassi Sergia_tenuiremis 0.01682182
## 28 Challengerosergia_hansjacobi Challengerosergia_talismani 1.00000000
## 29 Challengerosergia_hansjacobi Deosergestes_corniculum 1.00000000
## 30 Challengerosergia_hansjacobi Deosergestes_henseni 0.99999989
## 31 Challengerosergia_hansjacobi Eusergestes_arcticus 1.00000000
## 32 Challengerosergia_hansjacobi Gardinerosergia_splendens 1.00000000
## 33 Challengerosergia_hansjacobi Parasergestes_armatus 1.00000000
## 34 Challengerosergia_hansjacobi Parasergestes_vigilax 0.04921452
## 35 Challengerosergia_hansjacobi Phorcosergia_grandis 1.00000000
## 36 Challengerosergia_hansjacobi Robustasergia_robusta 1.00000000
## 37 Challengerosergia_hansjacobi Robustosergia_regalis 1.00000000
## 38 Challengerosergia_hansjacobi Sergestes_atlanticus 0.25203409
## 39 Challengerosergia_hansjacobi Sergia_tenuiremis 1.00000000
## 40 Challengerosergia_talismani Deosergestes_corniculum 0.99999999
## 41 Challengerosergia_talismani Deosergestes_henseni 0.99993757
## 42 Challengerosergia_talismani Eusergestes_arcticus 1.00000000
## 43 Challengerosergia_talismani Gardinerosergia_splendens 1.00000000
## 44 Challengerosergia_talismani Parasergestes_armatus 1.00000000
## 45 Challengerosergia_talismani Parasergestes_vigilax 0.02124209
## 46 Challengerosergia_talismani Phorcosergia_grandis 1.00000000
## 47 Challengerosergia_talismani Robustasergia_robusta 0.99999992
## 48 Challengerosergia_talismani Robustosergia_regalis 1.00000000
## 49 Challengerosergia_talismani Sergestes_atlanticus 0.11697344
## 50 Challengerosergia_talismani Sergia_tenuiremis 0.99999999
## 51 Deosergestes_corniculum Deosergestes_henseni 1.00000000
## 52 Deosergestes_corniculum Eusergestes_arcticus 0.99999999
## 53 Deosergestes_corniculum Gardinerosergia_splendens 1.00000000
## 54 Deosergestes_corniculum Parasergestes_armatus 1.00000000
## 55 Deosergestes_corniculum Parasergestes_vigilax 0.00236873
## 56 Deosergestes_corniculum Phorcosergia_grandis 1.00000000
## 57 Deosergestes_corniculum Robustasergia_robusta 1.00000000
## 58 Deosergestes_corniculum Robustosergia_regalis 1.00000000
## 59 Deosergestes_corniculum Sergestes_atlanticus 0.01350778
## 60 Deosergestes_corniculum Sergia_tenuiremis 1.00000000
## 61 Deosergestes_henseni Eusergestes_arcticus 0.99999959
## 62 Deosergestes_henseni Gardinerosergia_splendens 1.00000000
## 63 Deosergestes_henseni Parasergestes_armatus 0.99999978
## 64 Deosergestes_henseni Parasergestes_vigilax 0.00128436
## 65 Deosergestes_henseni Phorcosergia_grandis 1.00000000
## 66 Deosergestes_henseni Robustasergia_robusta 1.00000000
## 67 Deosergestes_henseni Robustosergia_regalis 1.00000000
## 68 Deosergestes_henseni Sergestes_atlanticus 0.00731827
## 69 Deosergestes_henseni Sergia_tenuiremis 1.00000000
## 70 Eusergestes_arcticus Gardinerosergia_splendens 0.99999997
## 71 Eusergestes_arcticus Parasergestes_armatus 1.00000000
## 72 Eusergestes_arcticus Parasergestes_vigilax 0.99996770
## 73 Eusergestes_arcticus Phorcosergia_grandis 1.00000000
## 74 Eusergestes_arcticus Robustasergia_robusta 0.99999993
## 75 Eusergestes_arcticus Robustosergia_regalis 1.00000000
## 76 Eusergestes_arcticus Sergestes_atlanticus 1.00000000
## 77 Eusergestes_arcticus Sergia_tenuiremis 0.99999911
## 78 Gardinerosergia_splendens Parasergestes_armatus 1.00000000
## 79 Gardinerosergia_splendens Parasergestes_vigilax 0.00266011
## 80 Gardinerosergia_splendens Phorcosergia_grandis 1.00000000
## 81 Gardinerosergia_splendens Robustasergia_robusta 1.00000000
## 82 Gardinerosergia_splendens Robustosergia_regalis 1.00000000
## 83 Gardinerosergia_splendens Sergestes_atlanticus 0.01499340
## 84 Gardinerosergia_splendens Sergia_tenuiremis 1.00000000
## 85 Parasergestes_armatus Parasergestes_vigilax 0.06386567
## 86 Parasergestes_armatus Phorcosergia_grandis 1.00000000
## 87 Parasergestes_armatus Robustasergia_robusta 1.00000000
## 88 Parasergestes_armatus Robustosergia_regalis 1.00000000
## 89 Parasergestes_armatus Sergestes_atlanticus 0.31517941
## 90 Parasergestes_armatus Sergia_tenuiremis 1.00000000
## 91 Parasergestes_vigilax Phorcosergia_grandis 0.00320670
## 92 Parasergestes_vigilax Robustasergia_robusta 0.00240837
## 93 Parasergestes_vigilax Robustosergia_regalis 0.00220257
## 94 Parasergestes_vigilax Sergestes_atlanticus 1.00000000
## 95 Parasergestes_vigilax Sergia_tenuiremis 0.00481337
## 96 Phorcosergia_grandis Robustasergia_robusta 1.00000000
## 97 Phorcosergia_grandis Robustosergia_regalis 1.00000000
## 98 Phorcosergia_grandis Sergestes_atlanticus 0.01851043
## 99 Phorcosergia_grandis Sergia_tenuiremis 1.00000000
## 100 Robustasergia_robusta Robustosergia_regalis 1.00000000
## 101 Robustasergia_robusta Sergestes_atlanticus 0.01358845
## 102 Robustasergia_robusta Sergia_tenuiremis 1.00000000
## 103 Robustosergia_regalis Sergestes_atlanticus 0.01271145
## 104 Robustosergia_regalis Sergia_tenuiremis 1.00000000
## 105 Sergestes_atlanticus Sergia_tenuiremis 0.02662028
## TestStat
## 1 7.666232e-02
## 2 1.070239e+01
## 3 1.269776e+01
## 4 1.710663e+01
## 5 1.836679e+01
## 6 1.514543e+00
## 7 1.687873e+01
## 8 1.006766e+01
## 9 8.341304e-01
## 10 1.659541e+01
## 11 1.702417e+01
## 12 1.732280e+01
## 13 8.875795e-04
## 14 1.526295e+01
## 15 9.695878e+00
## 16 1.293861e+01
## 17 1.810626e+01
## 18 2.051006e+01
## 19 1.218474e+00
## 20 1.764469e+01
## 21 8.943571e+00
## 22 1.268181e+00
## 23 2.023537e+01
## 24 1.726619e+01
## 25 2.072006e+01
## 26 7.817618e-02
## 27 1.423248e+01
## 28 2.850713e-03
## 29 1.459201e+00
## 30 2.159382e+00
## 31 7.165782e-01
## 32 1.474497e+00
## 33 3.154991e-03
## 34 1.218978e+01
## 35 7.281210e-01
## 36 1.670749e+00
## 37 1.215204e+00
## 38 8.958462e+00
## 39 1.652888e+00
## 40 1.971056e+00
## 41 2.909088e+00
## 42 7.908634e-01
## 43 1.901819e+00
## 44 1.381453e-02
## 45 1.378961e+01
## 46 1.203054e+00
## 47 2.130775e+00
## 48 1.877144e+00
## 49 1.051505e+01
## 50 1.917202e+00
## 51 1.091006e-01
## 52 1.922892e+00
## 53 9.383239e-03
## 54 1.528670e+00
## 55 1.795778e+01
## 56 4.919054e-01
## 57 3.949879e-02
## 58 9.043887e-02
## 59 1.464885e+01
## 60 1.486450e-01
## 61 2.285326e+00
## 62 4.454856e-02
## 63 2.224689e+00
## 64 1.912552e+01
## 65 1.138002e+00
## 66 1.181780e-02
## 67 4.397057e-01
## 68 1.581216e+01
## 69 1.724888e-02
## 70 2.036867e+00
## 71 6.530716e-01
## 72 2.808168e+00
## 73 1.445758e+00
## 74 2.123385e+00
## 75 1.716984e+00
## 76 1.501074e+00
## 77 2.362523e+00
## 78 1.546238e+00
## 79 1.773675e+01
## 80 5.367445e-01
## 81 9.239058e-03
## 82 1.427310e-01
## 83 1.445085e+01
## 84 8.641470e-02
## 85 1.169009e+01
## 86 7.995621e-01
## 87 1.740929e+00
## 88 1.285411e+00
## 89 8.475706e+00
## 90 1.723105e+00
## 91 1.738089e+01
## 92 1.792616e+01
## 93 1.809640e+01
## 94 6.837798e-01
## 95 1.660823e+01
## 96 7.055402e-01
## 97 2.447817e-01
## 98 1.405093e+01
## 99 7.801888e-01
## 100 2.376709e-01
## 101 1.463755e+01
## 102 4.529580e-02
## 103 1.476415e+01
## 104 3.827565e-01
## 105 1.336087e+01
##
## ------------------------------------------------------------
## Coefficients by group in variable "genus_species"
##
## Group: Allosergestes_pectinatus
## elevation slope
## estimate -2.812398 1.940053
## lower limit -3.491769 1.485663
## upper limit -2.133027 2.533418
##
## H0 : variables uncorrelated.
## R-squared : 0.8156794
## P-value : 9.643e-06
##
## Group: Allosergestes_sargassi
## elevation slope
## estimate -2.728530 1.844597
## lower limit -3.402995 1.419494
## upper limit -2.054065 2.397008
##
## H0 : variables uncorrelated.
## R-squared : 0.4037707
## P-value : 2.3852e-05
##
## Group: Challengerosergia_hansjacobi
## elevation slope
## estimate -1.528299 1.047036
## lower limit -1.913161 0.824145
## upper limit -1.143436 1.330207
##
## H0 : variables uncorrelated.
## R-squared : 0.4129122
## P-value : 3.3885e-06
##
## Group: Challengerosergia_talismani
## elevation slope
## estimate -1.578063 1.0389917
## lower limit -1.840235 0.8845412
## upper limit -1.315891 1.2204109
##
## H0 : variables uncorrelated.
## R-squared : 0.8120894
## P-value : 2.0471e-12
##
## Group: Deosergestes_corniculum
## elevation slope
## estimate -1.432046 0.8653057
## lower limit -1.737583 0.7020666
## upper limit -1.126509 1.0664999
##
## H0 : variables uncorrelated.
## R-squared : 0.8421211
## P-value : 8.1768e-08
##
## Group: Deosergestes_henseni
## elevation slope
## estimate -1.306994 0.8238705
## lower limit -1.593941 0.6637531
## upper limit -1.020047 1.0226131
##
## H0 : variables uncorrelated.
## R-squared : 0.3380975
## P-value : 1.6982e-06
##
## Group: Eusergestes_arcticus
## elevation slope
## estimate -2.1539322 1.3306114
## lower limit -3.8635100 0.6694485
## upper limit -0.4443543 2.6447543
##
## H0 : variables uncorrelated.
## R-squared : 0.8368056
## P-value : 0.029484
##
## Group: Gardinerosergia_splendens
## elevation slope
## estimate -1.2253950 0.8523703
## lower limit -1.5326713 0.6732695
## upper limit -0.9181186 1.0791149
##
## H0 : variables uncorrelated.
## R-squared : 0.4597308
## P-value : 1.1132e-06
##
## Group: Parasergestes_armatus
## elevation slope
## estimate -1.697577 1.0573117
## lower limit -2.095490 0.8217616
## upper limit -1.299663 1.3603801
##
## H0 : variables uncorrelated.
## R-squared : 0.5670657
## P-value : 1.5744e-06
##
## Group: Parasergestes_vigilax
## elevation slope
## estimate -3.306489 2.375713
## lower limit -4.455793 1.641298
## upper limit -2.157186 3.438750
##
## H0 : variables uncorrelated.
## R-squared : 0.5274941
## P-value : 0.00096145
##
## Group: Phorcosergia_grandis
## elevation slope
## estimate -1.470665 0.9366565
## lower limit -1.632022 0.8480626
## upper limit -1.309308 1.0345055
##
## H0 : variables uncorrelated.
## R-squared : 0.8798127
## P-value : < 2.22e-16
##
## Group: Robustasergia_robusta
## elevation slope
## estimate -1.2073093 0.8386026
## lower limit -1.5732188 0.6553974
## upper limit -0.8413998 1.0730198
##
## H0 : variables uncorrelated.
## R-squared : 0.6989428
## P-value : 6.7833e-07
##
## Group: Robustosergia_regalis
## elevation slope
## estimate -1.324440 0.8978007
## lower limit -1.531510 0.7807317
## upper limit -1.117369 1.0324239
##
## H0 : variables uncorrelated.
## R-squared : 0.8280327
## P-value : 2.4971e-15
##
## Group: Sergestes_atlanticus
## elevation slope
## estimate -2.839127 1.951572
## lower limit -3.727822 1.418514
## upper limit -1.950431 2.684945
##
## H0 : variables uncorrelated.
## R-squared : 0.7067609
## P-value : 8.6609e-05
##
## Group: Sergia_tenuiremis
## elevation slope
## estimate -1.2313110 0.8025570
## lower limit -1.6822103 0.5747311
## upper limit -0.7804117 1.1206941
##
## H0 : variables uncorrelated.
## R-squared : 0.2880646
## P-value : 0.0032339#save coefficients of fits as object
cc_species <- data.frame(coef(sma_species))This figure shows the static allometry of all the species with sufficient sampling (n = 15). It’s interactive, so you can click on different species in the legend to hide or show those points. That can be helpful for exploring patterns in a big mess of data. To me, it looks a bit like the smaller-bodied species have higher slopes? It looks like Sergestes atlanticus, Allosergestes_pectinatus, and Parasergestes vigilax have different (higher) slopes compared to the other sergestids, and the other species differ a bit in slopes but more so in intercepts. You can probably pull out which species are driving the difference by looking across the pairwise comparisons.
#reorder factor levels for figure legend (phylo order)
specimens_test$species_text <- factor(specimens_test$species_text,
levels = c("Deosergestes corniculum",
"Deosergestes henseni",
"Allosergestes pectinatus",
"Allosergestes sargassi",
"Sergestes atlanticus",
"Parasergestes vigilax",
"Parasergestes armatus",
"Eusergestes arcticus",
"Gardinerosergia splendens",
"Robustosergia regalis",
"Robustasergia robusta",
"Phorcosergia grandis",
"Sergia tenuiremis",
"Challengerosergia talismani",
"Challengerosergia hansjacobi"))
#make shape pallette
shapes.sp <- c("Deosergestes corniculum" = 21,
"Deosergestes henseni" = 22,
"Allosergestes pectinatus" = 23,
"Allosergestes sargassi" = 24,
"Sergestes atlanticus" = 25,
"Parasergestes vigilax" = 21,
"Parasergestes armatus" = 22,
"Eusergestes arcticus" = 23,
"Gardinerosergia splendens" = 24,
"Robustosergia regalis" = 25,
"Robustasergia robusta" = 21,
"Phorcosergia grandis" = 22,
"Sergia tenuiremis" = 23,
"Challengerosergia talismani" = 24,
"Challengerosergia hansjacobi" = 25)
cols.sp <- c("Deosergestes corniculum" = "#16ABA8",
"Deosergestes henseni" = "#1f78b4",
"Allosergestes pectinatus" = "#3acf75",
"Allosergestes sargassi" = "#33a02c",
"Sergestes atlanticus" = "#fb9a99",
"Parasergestes vigilax" = "#e31a1c",
"Parasergestes armatus" = "#fdbf6f",
"Eusergestes arcticus" = "#ff7f00",
"Gardinerosergia splendens" = "#cab2d6",
"Robustosergia regalis" = "#6a3d9a",
"Robustasergia robusta" = "gray31",
"Phorcosergia grandis" = "burlywood4",
"Sergia tenuiremis" = "black",
"Challengerosergia talismani" = "maroon1",
"Challengerosergia hansjacobi" = "orchid1")
#set s limits by species for plots
limits <- specimens_test %>%
group_by(genus_species) %>%
summarise(max = max(Body_Length, na.rm = TRUE),
min = min(Body_Length, na.rm = TRUE))
#make plot
plot_species <- ggplot(specimens_test, aes(x = Body_Length, y = Eye_Diameter, color = species_text, shape = species_text, fill = species_text)) +
geom_point(size = 2, alpha = 0.35) +
scale_shape_manual(values = shapes.sp, name = "Species") +
scale_color_manual(values = cols.sp, name = "Species") +
scale_fill_manual(values = cols.sp, name = "Species") +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.text = element_text(face = "italic")) +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_pectinatus")], xend = limits$max[which(limits$genus_species == "Allosergestes_pectinatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"])), color = cols.sp["Allosergestes pectinatus"]) + #Rana temporaria adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_sargassi")], xend = limits$max[which(limits$genus_species == "Allosergestes_sargassi")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"])), color = cols.sp["Allosergestes sargassi"]) + #Allosergestes_sargassi adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")], xend = limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"])), color = cols.sp["Challengerosergia hansjacobi"]) + #Challengerosergia_hansjacobi adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_talismani")], xend = limits$max[which(limits$genus_species == "Challengerosergia_talismani")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"])), color = cols.sp["Challengerosergia talismani"]) + #Challengerosergia_talismani adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_corniculum")], xend = limits$max[which(limits$genus_species == "Deosergestes_corniculum")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"])), color = cols.sp["Deosergestes corniculum"]) + #Deosergestes_corniculum adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_henseni")], xend = limits$max[which(limits$genus_species == "Deosergestes_henseni")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"])), color = cols.sp["Deosergestes henseni"]) + #Deosergestes henseni adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Gardinerosergia_splendens")], xend = limits$max[which(limits$genus_species == "Gardinerosergia_splendens")], y = 10^(log10(limits$min[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"])), color = cols.sp["Gardinerosergia splendens"]) + #Gardinerosergia splendens adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_armatus")], xend = limits$max[which(limits$genus_species == "Parasergestes_armatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"])), color = cols.sp["Parasergestes armatus"]) + #Parasergestes armatus adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_vigilax")], xend = limits$max[which(limits$genus_species == "Parasergestes_vigilax")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"])), color = cols.sp["Parasergestes vigilax"]) + #Parasergestes vigilax adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Phorcosergia_grandis")], xend = limits$max[which(limits$genus_species == "Phorcosergia_grandis")], y = 10^(log10(limits$min[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"])), color = cols.sp["Phorcosergia grandis"]) + #Phorcosergia grandis adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustasergia_robusta")], xend = limits$max[which(limits$genus_species == "Robustasergia_robusta")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustasergia_robusta")])*cc_species["Robustasergia_robusta", "slope"] + cc_species["Robustasergia_robusta", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustasergia_robusta")])*cc_species["Robustasergia_robusta", "slope"] + cc_species["Robustasergia_robusta", "elevation"])), color = cols.sp["Robustasergia robusta"]) + #Robustasergia robusta adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustosergia_regalis")], xend = limits$max[which(limits$genus_species == "Robustosergia_regalis")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"])), color = cols.sp["Robustosergia regalis"]) +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergestes_atlanticus")], xend = limits$max[which(limits$genus_species == "Sergestes_atlanticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"])), color = cols.sp["Sergestes atlanticus"]) +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergia_tenuiremis")], xend = limits$max[which(limits$genus_species == "Sergia_tenuiremis")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"])), color = cols.sp["Sergia tenuiremis"])+
geom_segment(aes(x = limits$min[which(limits$genus_species == "Eusergestes_arcticus")], xend = limits$max[which(limits$genus_species == "Eusergestes_arcticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"])), color = cols.sp["Eusergestes arcticus"])
#interactive plot
ggplotly(plot_species)#make plot with alpha = 1 for legend
plot_species_leg <- ggplot(specimens_test, aes(x = Body_Length, y = Eye_Diameter, color = species_text, shape = species_text, fill = species_text)) +
geom_point(size = 2, alpha = 1) +
scale_shape_manual(values = shapes.sp, name = "Species") +
scale_color_manual(values = cols.sp, name = "Species") +
scale_fill_manual(values = cols.sp, name = "Species") +
theme_bw()
# remove legend from previous plot
plot <- plot_species + theme(legend.position = "none")
# extract legend from this plot with alpha = 1
leg <- get_legend(plot_species_leg)
# construct full figure with transparent fig and non-transparent legend
plot_grid(plot, leg, ncol = 2, rel_widths = c(1, .4))
#export figure
#png("../Figures/dev-species.png", width = 12, height = 12, units = "in", res = 500)
pdf("../Figures/static-species.pdf", width = 9, height = 5)
plot_grid(plot, leg, ncol = 2, rel_widths = c(1, .4))
dev.off()Phylogenetic distribution of allometric slopes
Below, we plot species-specific slopes for developmental eye-body allometry onto the sergestid tree. There is potential to do more with these data; do we have predictions for what might drive species differences in developmental allometry? It looks like body size might be one factor, as the smaller-bodied species seem to have higher slopes. Whether species exhibit ontogenetic descent may be another factor to explore? We might predict that species that have ontogenetic descent and are going into dimmer habitats as they grow may invest in higher eye growth relative to body growth. Do all sergestids exhibit ontogenetic descent? So far, all I’ve done with these data in an evolutionary framework is test for phylogenetic signal. Let me know if you think of anything else cool to look at! Maybe we could test whether larger body sizes are correlated with lower slopes? Seemed from the plot of all the species that maybe the smaller bodied species had the higher slopes..
#Add slopes for developmental allometry into species dataframe ----
#add coefficient rownames to genus_species column
cc_species$genus_species <- rownames(cc_species)
#merge coefficients with species data by rowname
species_sma <- left_join(traits, cc_species, by = "genus_species") %>%
filter(genus_species != "Neosergestes_edwardsii")
# Prep data and phylogeny -----
#Make list of taxa to drop (in tree but not in dataset)
drops <- setdiff(tree$tip.label, as.character(species_sma$Binomial))
#Drop unwanted tips from phylogeny
tree.dev <- drop.tip(phy = tree, tip = drops)
# Make Genus_species tip labels to replace Binomial codes
sp.tips <- data.frame(old = species_sma$Binomial, new = species_sma$genus_species)
# replace tip labels in phylogeny
tree.dev.plot <- sub.taxa.label(tree.dev, sp.tips)
#change row names of species data
rownames(species_sma) <- species_sma$genus_species
#check that phylogeny and data match exactly
name.check(tree.dev.plot, species_sma)
#resort trait dataset to the order of tree tip labels
species_sma <- species_sma[tree.dev.plot$tip.label, ]
#make trait vector for slope
slope <- as.vector(species_sma$slope)
#add tip label names to vector
names(slope) <- species_sma$genus_species
# Phylogeny with slopes for developmental eye-body allometry----
plotTree.wBars(tree.dev.plot, slope,
scale = 0.5,
tip.labels = TRUE,
offset = 0.1,
ftype = "bi",
fsize = 1)
#export figure
#png("../Figures/phylo-slopes.png", width = 10, height = 7, units = "in", res = 500)
pdf("../Figures/phylo-slopes.pdf", width = 10, height = 7)
plotTree.wBars(tree.dev.plot, slope,
scale = 0.5,
tip.labels = TRUE,
offset = 0.1,
ftype = "bi",
fsize = 1)
dev.off()To test for phylogenetic signal, we fit Pagel’s lambda using a maximum likelihood approach in caper v.1.0.1.
#fit Pagel's lambda in caper
slopes <- data.frame(slope, names(slope))
mod <- pgls(slope ~ 1, comparative.data(tree.dev.plot, slopes,"names.slope."), lambda="ML")
#summary of model
summary(mod)##
## Call:
## pgls(formula = slope ~ 1, data = comparative.data(tree.dev.plot,
## slopes, "names.slope."), lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.97761 -0.02774 0.04743 0.33125 0.87950
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.770
## lower bound : 0.000, p = 0.20171
## upper bound : 1.000, p = 0.63856
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.30332 0.19267 6.7644 9.105e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5198 on 14 degrees of freedom
## Multiple R-squared: 0, Adjusted R-squared: 0
## F-statistic: NaN on 0 and 14 DF, p-value: NA#plot likelihood of lambda
plot(pgls.profile(mod, "lambda"))
The estimate for Pagel’s lambda is 0.77. However, this estimate was not significant, and neither the upper or lower bound of the confidence interval could be estimated. This could indicate that there is no phylogenetic signal in these data, but could also be the result of our sample size of species for comparative evolutionary models (n = 15). Simulations have previously found that estimating Pagel’s lambda for datasets of >20 species is problematic and often unreliable (see Freckleton et al., 2002). The likelihood profile of lambda shows that there is increased likelihood that it’s somewhere between 0.5 and 1.
Evolutionary eye-body allometry in Sergestidae
We next examine evolutionary (interspecific) eye-body allometry across sergestids using species means for eye diameter and body length. In calculating species means, we want to use adults only and also make sure that the data we use downstream for traits is relevant to the size class of individuals that we are taking species means from. The biggest concern for the Sergestidae is ontogenetic descent, as this is common across midwater animals and could confound relationships between eye size and depth. However, body size for ‘adult’ depth distribution data we used in this paper was cutoff between 4-11mm depending on species, and our smallest specimen here is 13 mm, so we are confident that species averages based on full species samples in our dataset are appropriate here.
Below, we use phylogenetic least-squares regression in caper to fit the allometric relationship between species means for eye size and body length across 16 species of sergestids.
# Prep data ------
#Calculate species means (3 outliers removed)
sp_means <- specimens_ed %>%
mutate_if(is.character, as.factor) %>%
group_by(genus_species) %>%
summarise(eye_av = mean(Eye_Diameter),
length_av = mean(Body_Length),
n = n()) %>%
ungroup()
#Merge with species trait data stored in dataframe 'species'
species <- traits %>%
left_join(sp_means, by = "genus_species") %>%
mutate(Organ = factor(Organ, levels = c("pesta","lensed", "unlensed", "none")))
#check that names match in dataframe and tree
name.check(phy = tree, data = species, data.names = species$Binomial)
#use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
species.comp <- comparative.data(data = species, phy = tree, names.col = "Binomial", vcv = TRUE, na.omit = FALSE, warn.dropped = TRUE)
#check for dropped tips or dropped species
species.comp$dropped$tips #phylogeny
species.comp$dropped$unmatched.rows #dataset
# PGLS model -------
#fit model for eye diameter ~ body length
pgls_evo <- pgls(log10(eye_av) ~ log10(length_av),
data = species.comp,
lambda = "ML", #uses Maximum Liklihood estimate of lambda
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_evo)
par(mfrow = c(1, 1))The standard model checks for the PGLS of interspecific eye diameter vs. body length are a bit wonky. The residuals look a little bimodal when they should be normal (though don’t have huge weird outliers). We’ll have to decide if we’re ok with this. Assuming we accept this, model results are below.
#print model output
summary(pgls_evo)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av), data = species.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.100464 -0.041912 0.009141 0.044272 0.080878
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.077147
## upper bound : 1.000, p = 1
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.288069 0.147956 -8.7057 5.051e-07 ***
## log10(length_av) 0.822852 0.095634 8.6042 5.810e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05855 on 14 degrees of freedom
## Multiple R-squared: 0.841, Adjusted R-squared: 0.8296
## F-statistic: 74.03 on 1 and 14 DF, p-value: 5.81e-07#save coefficients of fits as object
cc_pgls_evo <- data.frame(coef(pgls_evo))The PGLS regression shows that eyes scale with body length with a slope of 0.82, indicating that eyes have a negative evolutionary allometry (hypoallometric, slope < 1) with body length. This is what is most common among all vertebrates and some invertebrates that have been examined so far (inverts tend to vary a lot more than verts).
#Likelihood plot for Pagel's lambda
lambda.evo <- pgls.profile(pgls_evo, "lambda")
plot(lambda.evo, main = "Pagel's lambda",
sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")
The maximum likelihood estimate of lambda is 1.00 (high phylogenetic signal/Brownian Motion), so that’s what caper has used to create the variance-covariance matrix for the regression. However, small sample sizes (< 20-30 data points) have no/little power to estimate lambda (see Freckleton et al. 2002), so tests for phylogenetic signal in this dataset are essentially meaningless (note that the 95% confidence interval for lambda here is 0 to 1, all possible values, and the p-values are non-signficant). For this reason, if any of our other comparative analyses estimate a lambda of 0, we should repeat them again forcing lambda to 1, just so we can see what difference that makes in the model fits (as we can’t know what lambda actually is in this small dataset).
However, if we choose to keep this model, we can plot the model fit onto our data.
# Make plot
plot_evo <- ggplot(species, aes(x = Body_length, y = Eye_diameter, text = genus_species)) +
geom_point(size = 2, alpha = 0.6) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(data = cc_pgls_evo, aes(intercept = cc_pgls_evo[1,1], slope = cc_pgls_evo[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_pgls_evo[2,1], digits = 2), sep = " "), size = 3) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_pgls_evo[1,1], digits = 2), sep = " "), size = 3)
# Interactive plot
ggplotly(plot_evo)Bounded extremes of models
With sample sizes too small to generate a robust estimate of lambda, one option is to run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible parameter estimates are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.
Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_evo0 <- pgls(log10(eye_av) ~ log10(length_av),
data = species.comp,
lambda = 0, #set lambda to 0
bounds=list(lambda=c(0,0)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_evo0)
par(mfrow = c(1, 1))
#parameter estimates
summary(pgls_evo0)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av), data = species.comp,
## lambda = 0, param.CI = 0.95, bounds = list(lambda = c(0,
## 0)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.108877 -0.031123 -0.008501 0.040904 0.086975
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.431224 0.141272 -10.131 7.949e-08 ***
## log10(length_av) 0.928475 0.091351 10.164 7.636e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05883 on 14 degrees of freedom
## Multiple R-squared: 0.8807, Adjusted R-squared: 0.8721
## F-statistic: 103.3 on 1 and 14 DF, p-value: 7.636e-08#save coefficients of fits as object
cc_pgls_evo0 <- data.frame(coef(pgls_evo0))Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals). *Note that when we used a ML estimate of lambda, this was the value used, though the lambda estimate was not significant.
#fit model for eye diameter ~ body length (lambda = 1)
pgls_evo1 <- pgls(log10(eye_av) ~ log10(length_av),
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_evo1)
par(mfrow = c(1, 1))
#parameter estimates
summary(pgls_evo1)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av), data = species.comp,
## lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.100464 -0.041912 0.009141 0.044272 0.080878
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.288069 0.147956 -8.7057 5.051e-07 ***
## log10(length_av) 0.822852 0.095634 8.6042 5.810e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05855 on 14 degrees of freedom
## Multiple R-squared: 0.841, Adjusted R-squared: 0.8296
## F-statistic: 74.03 on 1 and 14 DF, p-value: 5.81e-07#save coefficients of fits as object
cc_pgls_evo1 <- data.frame(coef(pgls_evo1))We can plot both estimates of the relationship between eye size and body size onto our data to see the range of possible fits depending on the true value of lambda.
# Make plot
plot_evo_bounds <- ggplot(species, aes(x = Body_length, y = Eye_diameter, text = genus_species)) +
geom_point(size = 2, alpha = 0.6) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_abline(data = cc_pgls_evo0, aes(intercept = cc_pgls_evo0[1,1], slope = cc_pgls_evo0[2,1]), linetype = "solid") +
geom_abline(data = cc_pgls_evo1, aes(intercept = cc_pgls_evo1[1,1], slope = cc_pgls_evo1[2,1]), linetype = "dashed")
# Interactive plot
ggplotly(plot_evo_bounds)Eye diameter vs. body length across 16 species of sergesteds. Solid line indicates the PGLS fit when lambda = 0 and dashed line the fit when lambda = 1.
We can also plot the static allometry for each species onto the evolutionary allometry (across species means). Most of the species have very similar static/intraspecific slopes to the evolutionary/interspecific pattern, but those few small-bodied species are the strange ones here too.
# Make plot
plot_evo_bounds2 <- ggplot(species, aes(x = Body_length, y = Eye_diameter, text = genus_species)) +
geom_point(size = 2, alpha = 0.6) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_pectinatus")], xend = limits$max[which(limits$genus_species == "Allosergestes_pectinatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"])), color = "grey85") + #Allosergestes pectinatus adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_sargassi")], xend = limits$max[which(limits$genus_species == "Allosergestes_sargassi")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"])), color = "grey85") + #Allosergestes_sargassi adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")], xend = limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"])), color = "grey85") + #Challengerosergia_hansjacobi adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_talismani")], xend = limits$max[which(limits$genus_species == "Challengerosergia_talismani")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"])), color = "grey85") + #Challengerosergia_talismani adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_corniculum")], xend = limits$max[which(limits$genus_species == "Deosergestes_corniculum")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"])), color = "grey85") + #Deosergestes_corniculum adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_henseni")], xend = limits$max[which(limits$genus_species == "Deosergestes_henseni")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"])), color = "grey85") + #Deosergestes henseni adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Gardinerosergia_splendens")], xend = limits$max[which(limits$genus_species == "Gardinerosergia_splendens")], y = 10^(log10(limits$min[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"])), color = "grey85") + #Gardinerosergia splendens adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_armatus")], xend = limits$max[which(limits$genus_species == "Parasergestes_armatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"])), color = "grey85") + #Parasergestes armatus adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_vigilax")], xend = limits$max[which(limits$genus_species == "Parasergestes_vigilax")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"])), color = "grey85") + #Parasergestes vigilax adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Phorcosergia_grandis")], xend = limits$max[which(limits$genus_species == "Phorcosergia_grandis")], y = 10^(log10(limits$min[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"])), color = "grey85") + #Phorcosergia grandis adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustasergia_robusta")], xend = limits$max[which(limits$genus_species == "Robustasergia_robusta")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustasergia_robusta")])*cc_species["Robustasergia_robusta", "slope"] + cc_species["Robustasergia_robusta", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustasergia_robusta")])*cc_species["Robustasergia_robusta", "slope"] + cc_species["Robustasergia_robusta", "elevation"])), color = "grey85") + #Robustasergia robusta adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustosergia_regalis")], xend = limits$max[which(limits$genus_species == "Robustosergia_regalis")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"])), color = "grey85") +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergestes_atlanticus")], xend = limits$max[which(limits$genus_species == "Sergestes_atlanticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"])), color = "grey85") +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergia_tenuiremis")], xend = limits$max[which(limits$genus_species == "Sergia_tenuiremis")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"])), color = "grey85") +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Eusergestes_arcticus")], xend = limits$max[which(limits$genus_species == "Eusergestes_arcticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"])), color = "grey85") +
geom_abline(data = cc_pgls_evo0, aes(intercept = cc_pgls_evo0[1,1], slope = cc_pgls_evo0[2,1]), linetype = "solid") +
geom_abline(data = cc_pgls_evo1, aes(intercept = cc_pgls_evo1[1,1], slope = cc_pgls_evo1[2,1]), linetype = "dashed")
# Interactive plot
plot_evo_bounds2
#export evolutionary allometry figures to pdf
#witout static allom
pdf("../Figures/evo_allom.pdf", width = 4, height = 2.5)
plot_evo_bounds
dev.off()## quartz_off_screen
## 2#with static allom
pdf("../Figures/evo_allom-static.pdf", width = 4, height = 2.5)
plot_evo_bounds2
dev.off()## quartz_off_screen
## 2Correlates of eye size in deep-sea sergestid shrimp
Next, we examine whether variation in eye size is correlated with ecological variables related to light regime and vision: depth and bioluminescent organs.
Phylogenetic distribution of eye size and ecology
Here, we examine how absolute and relative eye size vary across the sergestid phylogeny, and the corresponding distribution ecological variables. For visualizing relative eye size, we use the residuals of the PGLS fit for eye-body allometry, as these show eye size corrected for both body size and evolutionary relationships. (Here I have used the model where lambda = 1, as that was the ML estimate despite being non-significant.)
#Add residuals from PGLS fits to dataframe ----
#extract pgls residuals
pglsres <- residuals(pgls_evo)
#name residuals
colnames(pglsres) <- "pglsres"
#add rownames to species dataframe
rownames(species) <- species$Binomial
#merge residuals with species data by rowname
species.phy <- merge(species, pglsres, by = "row.names")
#make column converted to observed eye size as factor of PGLS fit
species.phy$eyefactor <- 10^species.phy$pglsres
# Prep phylogeny -----
# Make Genus_species tip labels to replace Binomial codes
sp.tips <- data.frame(old = species.phy$Binomial, new = species.phy$genus_species)
# replace tip labels in phylogeny
tree.plot <- sub.taxa.label(tree, sp.tips)
# Prep data ----
#change row names of species data
rownames(species.phy) <- species.phy$genus_species
#check that phylogeny and data match exactly
name.check(tree.plot, species.phy)
#resort trait dataset to the order of tree tip labels
species.phy <- species.phy[tree.plot$tip.label, ]
#make trait vector for absolute eye size
aveye <- as.vector(species.phy$eye_av)
names(aveye) <- species.phy$genus_species
#make trait vector for eye investment (PGLS residuals)
releye <- as.vector(species.phy$pglsres) #residuals of pgls
names(releye) <- species.phy$genus_species
# Phylogeny with eye size and investment -------
#color vector for photophores
col_ph <- c("lensed" = "palegreen",
"none" = "honeydew",
"pesta" = "cadetblue1",
"unlensed" = "orchid1")
#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye,
scale = 1.5,
tip.labels = TRUE,
col = unname(col_ph[as.factor(species.phy$Organ)]),
offset = 1,
ftype = "bi",
fsize = 0.7)
#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
pch = 19, #shape of labels
col = "mediumpurple",
offset = 0)
#add legend for photophores
legend(x = "bottomleft", legend = c("lensed", "none", "pesta", "unlensed"), pch = 22, pt.cex= 2, pt.bg = col_ph, cex = 1, bty = "n", horiz = F)
Phylogeny of the 15 sertestids sampled depicting absolute eye diameter (purple circles scaled by eye diameter), relative eye investment (bars scaled by the residuals from the PGLS of eye diameter vs. body length, lambda = 1), and photophore complexity (color of bars).
#Export figure------
pdf("../Figures/phylo-photo1.pdf", width = 6, height = 4)
#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye,
scale = 1.5,
tip.labels = TRUE,
col = unname(col_ph[species.phy$Organ]),
offset = 1,
ftype = "bi",
fsize = 0.7)
#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
pch = 19, #shape of labels
col = "mediumpurple",
offset = 0)
#add legend for photophores
legend(x = "bottomleft", legend = c("lensed", "none", "pesta", "unlensed"), pch = 22, pt.cex= 2, pt.bg = col_ph, cex = 1, bty = "n", horiz = F)
dev.off()We can also visualize day and night depth distributions.
# tree showing day and night mean depths across species
#fix node order upset by ladderization
#filter out internal nodes from second column of edge matrix
is_tip <- tree.plot$edge[,2] <= length(tree.plot$tip.label)
ordered_tips <- tree.plot$edge[is_tip, 2]
#extract tips in right order
#tree.plot$tip.label[ordered_tips]
#resort trait dataset to the order of tree tip labels
species.phy <- species.phy[tree.plot$tip.label[ordered_tips], ]
#make trait vector for day depth av
day <- as.vector(species.phy$Day_Avg)
names(day) <- species.phy$genus_species
#make trait vector for night depth av
night <- as.vector(species.phy$Night_Avg)
names(night) <- species.phy$genus_species
#make trait vector for max reported depth
max <- as.vector(species.phy$Max_reported)
names(max) <- species.phy$genus_species
#Plot tree with depth ranges
plot(tree.plot,
show.tip.label= TRUE,
font=3,
cex=0.9,
label.offset = 0.8,
edge.width=2)
#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
pch = 19, #shape of labels
col = "mediumpurple",
offset = 0)
#add gray bar for difference in day/night depths (migration distance)
segments(1.05 + (day*0.0003), 1:length(day), 1.05 + (night*0.0003), 1:length(night), col="gray", lwd=8)
#add yellow point for day depth
segments(1.05 + (day*0.0003), 1:length(day), 1.05 + (day*0.0003), 1:length(day), col="orange", lwd=10)
#add blue point for night depth
segments(1.05 + (night*0.0003), 1:length(night), 1.05 + (night*0.0003), 1:length(night), col="royalblue2", lwd=10)
#add black point for max depth
segments(1.05 + (max*0.0003), 1:length(max), 1.05 + (max*0.0003), 1:length(max), col="black", lwd=7)
# export figure
pdf("../Figures/phylo-depth.pdf", width = 8, height = 5.5)
#Plot tree with depth ranges
plot(tree.plot,
show.tip.label= TRUE,
font=3,
cex=0.9,
label.offset = 0.8,
edge.width=2)
#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
pch = 19, #shape of labels
col = "mediumpurple",
offset = 0)
#add gray bar for difference in day/night depths (migration distance)
segments(1.05 + (day*0.0003), 1:length(day), 1.05 + (night*0.0003), 1:length(night), col="gray", lwd=8)
#add yellow point for day depth
segments(1.05 + (day*0.0003), 1:length(day), 1.05 + (day*0.0003), 1:length(day), col="orange", lwd=10)
#add blue point for night depth
segments(1.05 + (night*0.0003), 1:length(night), 1.05 + (night*0.0003), 1:length(night), col="royalblue2", lwd=10)
#add black point for max depth
segments(1.05 + (max*0.0003), 1:length(max), 1.05 + (max*0.0003), 1:length(max), col="black", lwd=7)
dev.off()The above plot shows absolute eye size (purple dots scaled by eye diameter) and the average daytime (orange) and nighttime (blue) depths for each species in the phylogeny, as well as the maximum reported depth (black). Gray bars indicate vertical distance traveled during diel vertical migrations. Note that this plot is missing a scale bar, easier to add in Illustrator, but minimum depth at night (left/blue points) is 133m and maximum depth during the day (right/orange points) is 1750m. The maximum depth of occurrance (black points) across species is 2300m.
Photophore complexity
Does eye size correlate with photophore complexity?
Prediction: Eye size may be largest in animals with unlensed photophores, followed by lensed and organ of pesta types, as organs that are visually more difficult to resolve will require greater sensitivity to light (note: should this say longer focal lengths or something like that for higer acuity, not sensitivity?).
Preliminary plots (photophores)
Before incorporating phylogeny, we can inspect the distribution of our data with boxplots for absolute and relative eye size by photophore type.
#Boxplots for absolute eye size across photophores -----
plot_abs_photo <- ggplot(data = species.phy, aes(x = reorder(Organ, eye_av, FUN = mean), y = eye_av, text = genus_species)) +
geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
geom_jitter(shape = 19, size = 1.3, alpha = 0.5, position = position_jitter(0.15)) +
theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black")) + #controls background
xlab("Photophore complexity") +
ylab("Eye diameter (mm)") +
ggtitle("Absolute eye size")
ggplotly(plot_abs_photo)#Boxplots for relative eye investment across photophores -----
plot_rel_photo <- ggplot(data = species.phy, aes(x = reorder(Organ, pglsres, FUN = mean), y = pglsres, text = genus_species)) +
geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
geom_jitter(shape = 19, size = 1.3, alpha = 0.5, position = position_jitter(0.15)) +
theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black")) + #controls background
xlab("Photophore complexity") +
ylab("Eye ~ body residuals") +
ggtitle("Relative eye investment")
ggplotly(plot_rel_photo) It looks like species with organs of pesta tend to have smaller eye diameters and smaller relative eye sizes than species with other types of photophores (note: we may want to combine the non-pesta organs together for power as sample size within groups is really small). We next test this in an evolutionary framework.
Test absolute eye size across photophores (PGLS)
This model is for eye diameter vs. photophore type. Lambda is estimated by maximum likelihood.
#run PGLS model using the ML estimate of lambda
pgls_abs_photo <- pgls(eye_av ~ Organ,
data = species.comp,
lambda = "ML",
#uses Maximum Liklihood estimate of lambda
param.CI = 0.95)
#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_abs_photo) #plot the linear model
par(mfrow = c(1,1)) #set plotting window back to one plot
#main effects
anova(pgls_abs_photo)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Organ 3 0.71605 0.238683 2.8559 0.08159 .
## Residuals 12 1.00292 0.083576
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#coefficients
summary(pgls_abs_photo)##
## Call:
## pgls(formula = eye_av ~ Organ, data = species.comp, lambda = "ML",
## param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.55499 -0.18875 -0.03819 0.05956 0.27422
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.21126
## upper bound : 1.000, p = 1
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.87918 0.12205 7.2035 1.081e-05 ***
## Organlensed 0.21111 0.26870 0.7857 0.44729
## Organunlensed 0.56885 0.20812 2.7332 0.01816 *
## Organnone 0.32053 0.28563 1.1222 0.28373
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2891 on 12 degrees of freedom
## Multiple R-squared: 0.4166, Adjusted R-squared: 0.2707
## F-statistic: 2.856 on 3 and 12 DF, p-value: 0.08159#Likelihood plot for Pagel's lambda
plot(pgls.profile(pgls_abs_photo, "lambda"),
main = "Pagel's lambda",
sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")
This test shows no significant effect of photophore type (pesta, lensed, unlensed, none) on mean absolute eye size across species. The ML estimate of lambda in this model is 1, but this is not significant. So, next we fix lambda at its bounds (0 and 1) to see the effects on our results.
Test when lambda = 0
Absolute eye size vs. photophore type
#fit model for eye diameter ~ body length (lambda = 0)
pgls_abs_photo0 <- pgls(eye_av ~ Organ,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_abs_photo0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_abs_photo0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Organ 3 1.24698 0.41566 5.4302 0.01361 *
## Residuals 12 0.91855 0.07655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_abs_photo0)##
## Call:
## pgls(formula = eye_av ~ Organ, data = species.comp, lambda = 1e-05,
## param.CI = 0.95, bounds = list(lambda = c(1e-05, 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.49463 -0.23850 -0.09265 0.12102 0.48200
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.829177 0.092223 8.9910 1.116e-06 ***
## Organlensed 0.319005 0.216283 1.4749 0.165975
## Organunlensed 0.653861 0.166258 3.9328 0.001989 **
## Organnone 0.420822 0.291635 1.4430 0.174619
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2767 on 12 degrees of freedom
## Multiple R-squared: 0.5758, Adjusted R-squared: 0.4698
## F-statistic: 5.43 on 3 and 12 DF, p-value: 0.01361When lambda is set to 0, there is a significant effect of organ type on absolute eye size.
Test when lambda = 1
Absolute eye size vs. photophore type.
#fit model
pgls_abs_photo1 <- pgls(eye_av ~ Organ,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_abs_photo1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_abs_photo1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Organ 3 0.71605 0.238683 2.8559 0.08159 .
## Residuals 12 1.00292 0.083576
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_abs_photo1)##
## Call:
## pgls(formula = eye_av ~ Organ, data = species.comp, lambda = 1,
## param.CI = 0.95, bounds = list(lambda = c(1, 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.55499 -0.18875 -0.03819 0.05956 0.27422
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.87918 0.12205 7.2035 1.081e-05 ***
## Organlensed 0.21111 0.26870 0.7857 0.44729
## Organunlensed 0.56885 0.20812 2.7332 0.01816 *
## Organnone 0.32053 0.28563 1.1222 0.28373
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2891 on 12 degrees of freedom
## Multiple R-squared: 0.4166, Adjusted R-squared: 0.2707
## F-statistic: 2.856 on 3 and 12 DF, p-value: 0.08159When lambda is set to 1, there is no significant effect of photophore type (pesta, lensed, unlensed, none) on mean absolute eye diameter across species.
Note here that though our overall test for significant effect of photophore type (the anova) is not significant (p = 0.08), the contrasts show one significant difference between the unlensed and the pesta states. The contrasts here are more post-hoc, so we shouldn’t really consider them given that the full model was not significant. Full model wins out. For the paper, I’d only report these contrasts for models with a significant effect of the ecological trait being tested.
Test relative eye size across photophores (PGLS)
Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with photophore complexity as a covariate.
#run PGLS model using the Maximum Liklihood estimate of lambda
pgls_photo <- pgls(log10(eye_av) ~ log10(length_av) + Organ,
data = species.comp,
lambda = "ML",
#uses Maximum Liklihood estimate of lambda
param.CI = 0.95)
#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_photo) #plot the linear model
par(mfrow = c(1,1)) #set plotting window back to one plotThe model checks look pretty good here.
Next we can take a look at the overall significance of the main effects and sequential sums of squares. Here we can see that body length has a significant effect on eye size (expected, we already knew that), but we also find that photophore complexity/type has a significant effect on eye size when the effect of body size is accounted for (p = 0.01).
#anova
anova(pgls_photo)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.35757 0.35757 213.2573 1.511e-08 ***
## Organ 3 0.03002 0.01001 5.9672 0.01144 *
## Residuals 11 0.01844 0.00168
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Finally, we look at the PGLS coefficients as well as significance with respect to contrasts. Note that for a categorical variable, R will treat whichever factor comes first as the “control” and will test other factors for significant differences relative to that “control”. Here, the comparison state (“control”) for photophore complexity is set to “pesta” (estimate listed under (Intercept) - that’s why you don’t see that state as a comparison line) and the other states show the difference in the estimated intercept for that state compared to the pesta state, and whether that difference is significant. You can see that “unlensed” and “lensed” both show significant differences from “pesta”, and the intercept for both is higher, indicating larger mean relative eye sizes in those groups.
#summary
summary(pgls_photo)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Organ, data = species.comp,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.083912 -0.026000 -0.007001 0.014353 0.041549
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.099077
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.229881 0.114139 -10.7753 3.489e-07 ***
## log10(length_av) 0.769679 0.077417 9.9420 7.835e-07 ***
## Organlensed 0.099594 0.032644 3.0509 0.011029 *
## Organunlensed 0.108865 0.029393 3.7038 0.003478 **
## Organnone 0.051706 0.045672 1.1321 0.281667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04095 on 11 degrees of freedom
## Multiple R-squared: 0.9546, Adjusted R-squared: 0.9381
## F-statistic: 57.79 on 4 and 11 DF, p-value: 2.576e-07The output also includes the maximum-likelihood estimation of lambda, which represents the phylogenetic signal for PGLS model residuals. Here, the estimate for lambda is 0.00, which means that there is no phylogenetic signal. This doesn’t mean we shouldn’t include phylogeny in the model (no reason not to), but it does mean that it’s the same result we would get if we assessed this without including phylogeny. This also doesn’t mean that there is not phylogenetic signal in the individual traits that we put into the model; this is only lambda for the model residuals, not the traits themselves. However, as mentioned above, when you have <20 species the power to estimate lambda is not there, so we may want to consider running these models twice, one forcing lambda to 1 and another forcing lambda to 0, so that we can see how that affects our findings. That could help us in our interpretation (because photophore type is tightly matched to clades, my expectation is that if we force lambda to 0.5 or 1, the relationship may not remain significant.)
#Likelihood plot for Pagel's lambda
plot(pgls.profile(pgls_photo, "lambda"),
main = "Pagel's lambda",
sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")
With sample sizes too small to generate a robust estimate of lambda, one option is to run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible parameter estimates are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.
Test when lambda = 0
Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals). *Note that when we used a ML estimate of lambda, this was the value used, though the lambda estimate was not significant.
#fit model for eye diameter ~ body length (lambda = 0)
pgls_photo0 <- pgls(log10(eye_av) ~ log10(length_av) + Organ,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo0)
par(mfrow = c(1, 1))
#parameter estimates
summary(pgls_photo0)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Organ, data = species.comp,
## lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.083912 -0.026001 -0.007001 0.014353 0.041549
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.229880 0.114139 -10.7753 3.489e-07 ***
## log10(length_av) 0.769678 0.077417 9.9420 7.835e-07 ***
## Organlensed 0.099594 0.032644 3.0509 0.011029 *
## Organunlensed 0.108865 0.029393 3.7038 0.003478 **
## Organnone 0.051706 0.045673 1.1321 0.281667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04095 on 11 degrees of freedom
## Multiple R-squared: 0.9546, Adjusted R-squared: 0.9381
## F-statistic: 57.79 on 4 and 11 DF, p-value: 2.576e-07anova(pgls_photo0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.35757 0.35757 213.2555 1.511e-08 ***
## Organ 3 0.03002 0.01001 5.9671 0.01144 *
## Residuals 11 0.01844 0.00168
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#save coefficients of fits as object
cc_pgls_photo0 <- data.frame(coef(pgls_photo0))Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 1)
pgls_photo1 <- pgls(log10(eye_av) ~ log10(length_av) + Organ,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo1)
par(mfrow = c(1, 1))
#model outputs
anova(pgls_photo1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.253779 0.253779 106.0654 5.504e-07 ***
## Organ 3 0.021672 0.007224 3.0193 0.07574 .
## Residuals 11 0.026319 0.002393
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(pgls_photo1)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Organ, data = species.comp,
## lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.067903 -0.023288 -0.004012 0.039617 0.064583
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.185672 0.135069 -8.7783 2.674e-06 ***
## log10(length_av) 0.740092 0.089218 8.2953 4.620e-06 ***
## Organlensed 0.102270 0.045518 2.2468 0.04615 *
## Organunlensed 0.107974 0.038334 2.8167 0.01677 *
## Organnone 0.066800 0.049589 1.3471 0.20505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04891 on 11 degrees of freedom
## Multiple R-squared: 0.9128, Adjusted R-squared: 0.8811
## F-statistic: 28.78 on 4 and 11 DF, p-value: 8.972e-06Photophores have a significant effect on eye size when lambda is set to 0 (no phylogenetic signal in model residuals). However, we lose that significance when lambda is set to 1 (high phylogenetic signal in model residuals). While not significant here, it is close (p = 0.08). We can’t be totally confident in our finding, as it depends on the true value of lambda (which we can’t estimate with our sample size of species). However, since all patterns we see seem to be driven by organ of pesta, maybe we could try combining categories into organ of pesta or not, for reduced states and more power? I have done that next.
Photophore complexity (pesta vs. not - reduced states)
Here, I’ve reduced the number of states to the following: bioluminescent with organ of pesta present (“pesta”), bioluminescent with organ of pesta not present (“nonpesta”), bioluminescence not present (“none”). Then I’ve re-run the models with lambda set to 0 and 1 to see the binned range of potential findings depending on the true value of lambda.
#recategorize photophore bins
species2 <- species %>%
mutate(photo2 = recode(Organ, "lensed" = "nonpesta", "unlensed" = "nonpesta")) %>%
mutate(photo2 = factor(photo2, levels = c("pesta","nonpesta","none")))
#use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
species.comp2 <- comparative.data(data = species2, phy = tree, names.col = "Binomial", vcv = TRUE, na.omit = FALSE, warn.dropped = TRUE)
#check for dropped tips or dropped species
species.comp2$dropped$tips #phylogeny
species.comp2$dropped$unmatched.rows #datasetData visualization
First we can look at the distribution of data again in these new states with a phylogeny visualization and with boxplots.
# Phylogeny with eye size and investment -------
#color palette for photophore groups
col_photo2 <- c("pesta" = "#7b3294",
"lensed" = "#008837",
"unlensed" = "#a6dba0",
"none" = "#A6ACAF")
#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye,
scale = 1.5,
tip.labels = TRUE,
col = "gray",
offset = 1,
ftype = "bi",
fsize = 0.7)
#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
pch = 19, #shape of labels
col = unname(col_photo2[species2$Organ]),
offset = 0)
#add legend for photophores
legend(x = "bottomleft", legend = c("pesta", "lensed", "unlensed", "none"), pch = 22, pt.cex= 2, pt.bg = col_photo2, cex = 1, bty = "n", horiz = F)
Phylogeny of the 15 sertestids sampled depicting absolute eye diameter (circles, colored by photophore state) and relative eye investment (bars scaled by the residuals from the PGLS of eye diameter vs. body length, lambda = 1).
#Export figure------
pdf("../Figures/phylo-photo2.pdf", width = 6, height = 4)
#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye,
scale = 1.5,
tip.labels = TRUE,
col = "gray",
offset = 1,
ftype = "bi",
fsize = 0.7)
#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
pch = 19, #shape of labels
col = unname(col_photo2[species2$Organ]),
offset = 0)
#add legend for photophores
legend(x = "bottomleft", legend = c("pesta", "lensed", "unlensed", "none"), pch = 22, pt.cex= 2, pt.bg = col_photo2, cex = 1, bty = "n", horiz = F)
dev.off()sh_photo2 <- c("lensed" = 17,
"none" = 18,
"pesta" = 19,
"unlensed" = 15)
#Boxplots for absolute eye size across photophores -----
plot_abs_photo2 <- ggplot(data = species2, aes(x = factor(photo2, levels=c("pesta", "nonpesta","none")), y = eye_av)) +
geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
geom_jitter(aes(x = photo2, y = eye_av, color = Organ, shape = Organ), size = 3, alpha = 0.9, position = position_jitter(0.15)) +
theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black"), legend.key = element_blank()) + #controls background
scale_color_manual(values = col_photo2, name = "Photophore type",
breaks = c("pesta", "lensed", "unlensed", "none")) +
scale_shape_manual(values = sh_photo2, name = "Photophore type",
breaks = c("pesta", "lensed", "unlensed", "none")) +
xlab("Photophore complexity") +
ylab("Eye diameter (mm)") +
ggtitle("Absolute eye size")
plot_abs_photo2
#Boxplots for relative eye size across photophores -----
plot_rel_photo2 <- ggplot(data = species2, aes(x = factor(photo2, levels=c("pesta", "nonpesta","none")), y = eye_av/length_av)) +
geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
geom_jitter(aes(x = photo2, y = eye_av/length_av, color = Organ, shape = Organ), size = 3, alpha = 0.9, position = position_jitter(0.15)) +
theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black"), legend.key = element_blank()) + #controls background
scale_color_manual(values = col_photo2, name = "Photophore type",
breaks = c("pesta", "lensed", "unlensed", "none")) +
scale_shape_manual(values = sh_photo2, name = "Photophore type",
breaks = c("pesta", "lensed", "unlensed", "none")) +
xlab("Photophore complexity") +
ylab("Relative eye size (eye/body)") +
ggtitle("Relative eye size")
plot_rel_photo2
#Boxplots for residual eye investment across photophores -----
plot_inv_photo2 <- ggplot(data = species2, aes(x = factor(photo2, levels=c("pesta", "nonpesta","none")), y = pglsres)) +
geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
geom_jitter(aes(x = photo2, y = pglsres, color = Organ, shape = Organ), size = 3, alpha = 0.9, position = position_jitter(0.15)) +
theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black"), legend.key = element_blank()) + #controls background
scale_color_manual(values = col_photo2, name = "Photophore type",
breaks = c("pesta", "lensed", "unlensed", "none")) +
scale_shape_manual(values = sh_photo2, name = "Photophore type",
breaks = c("pesta", "lensed", "unlensed", "none")) +
xlab("Photophore complexity") +
ylab("Eye ~ body residuals") +
ggtitle("Residual eye investment")
plot_inv_photo2
# Make boxplot figure to export ----
#name panels
fig.a <- plot_abs_photo2 + theme(legend.position = "none") + ggtitle("") + xlab("")
fig.b <- plot_rel_photo2 + theme(legend.position = "none") + ggtitle("") + xlab("")
fig.c <- plot_inv_photo2 + theme(legend.position = "none") + ggtitle("") + xlab("")
#arrange plots in panels
plots <- plot_grid(fig.a, fig.b, fig.c, #list of plots to arrange in grid
align = 'vh', #align horizontally and vertically
axis = 'lb',
labels = c("A", "B", "C"), #panel labels for figure
hjust = -1, #adjustment for panel labels
nrow = 1) #number of rows in grids
# extract legend from Rana temporaria figure
leg <- get_legend(fig.a + theme(legend.position="right"))
#export figure
pdf("../Figures/boxplots-photo.pdf", width = 10, height = 4)
plot_grid(plots, leg, nrow = 1, rel_widths = c(1, .2))
dev.off()## quartz_off_screen
## 2Test absolute eye size vs. organ type (pesta, nonpesta, none)
This is a model of absolute eye diameter vs. organ type (PGLS). Again, we will first run a model that estimates lambda by ML.
#fit model
pgls_photo_bin_abs <- pgls(eye_av ~ photo2,
data = species.comp2,
lambda = "ML",
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo_bin_abs)
par(mfrow = c(1, 1))
#main effects
anova(pgls_photo_bin_abs)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## photo2 2 0.50453 0.252266 2.7004 0.1045
## Residuals 13 1.21443 0.093418#parameter estimates
summary(pgls_photo_bin_abs)##
## Call:
## pgls(formula = eye_av ~ photo2, data = species.comp2, lambda = "ML",
## param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.56167 -0.26508 -0.06907 0.06838 0.29391
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.33913
## upper bound : 1.000, p = 1
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.87918 0.12904 6.8135 1.237e-05 ***
## photo2nonpesta 0.49215 0.21405 2.2992 0.03871 *
## photo2none 0.33982 0.30171 1.1263 0.28039
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3056 on 13 degrees of freedom
## Multiple R-squared: 0.2935, Adjusted R-squared: 0.1848
## F-statistic: 2.7 on 2 and 13 DF, p-value: 0.1045When lambda is fitted by ML, the estimate for lambda is 1; however, this lambda estimate is not significant. Again, we can run two models bounding the extreme possibilities of lambda to see the range of possible results.
Test with lambda = 0
Absolute eye size vs. organ type (pesta, non-pesta, none)
#fit model for eye diameter ~ body length (lambda = 0)
pgls_photo0_bin_abs <- pgls(eye_av ~ photo2,
data = species.comp2,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo0_bin_abs)
par(mfrow = c(1, 1))
#main effects
anova(pgls_photo0_bin_abs)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## photo2 2 1.0975 0.54874 6.6791 0.01011 *
## Residuals 13 1.0680 0.08216
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_photo0_bin_abs)##
## Call:
## pgls(formula = eye_av ~ photo2, data = species.comp2, lambda = 1e-05,
## param.CI = 0.95, bounds = list(lambda = c(1e-05, 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.57427 -0.26319 -0.03726 0.13291 0.39016
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.829177 0.095544 8.6785 9.072e-07 ***
## photo2nonpesta 0.542243 0.151069 3.5894 0.003298 **
## photo2none 0.420823 0.302137 1.3928 0.187036
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2866 on 13 degrees of freedom
## Multiple R-squared: 0.5068, Adjusted R-squared: 0.4309
## F-statistic: 6.679 on 2 and 13 DF, p-value: 0.01011When lambda = 0, there is a significant effect of photophore type on absolute eye size, and the contrasts show that species with organs of pesta have significantly smaller mean relative eye sizes than those with other types of organ.
Test with lambda = 1
#fit model for eye diameter ~ body length (lambda = 1)
pgls_photo1_bin_abs <- pgls(eye_av ~ photo2,
data = species.comp2,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo1_bin_abs)
par(mfrow = c(1, 1))
#main effects
anova(pgls_photo1_bin_abs)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## photo2 2 0.50453 0.252266 2.7004 0.1045
## Residuals 13 1.21443 0.093418#parameter estimates
summary(pgls_photo1_bin_abs)##
## Call:
## pgls(formula = eye_av ~ photo2, data = species.comp2, lambda = 1,
## param.CI = 0.95, bounds = list(lambda = c(1, 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.56167 -0.26508 -0.06907 0.06838 0.29391
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.87918 0.12904 6.8135 1.237e-05 ***
## photo2nonpesta 0.49215 0.21405 2.2992 0.03871 *
## photo2none 0.33982 0.30171 1.1263 0.28039
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3056 on 13 degrees of freedom
## Multiple R-squared: 0.2935, Adjusted R-squared: 0.1848
## F-statistic: 2.7 on 2 and 13 DF, p-value: 0.1045When lambda = 1, there is no significant effect of photophore type on absolute eye size (this kind of makes sense, as all the species with the organs of pesta are one clade with one origin).
Test relative eye size vs. organ type
Here is a model of log10(eye diameter) vs. log10(body length) + photophore type (pesta, nonpesta, none), with lambda fitted by ML.
#fit model
pgls_photo_bin <- pgls(log10(eye_av) ~ log10(length_av) + photo2,
data = species.comp2,
lambda = "ML",
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo_bin)
par(mfrow = c(1, 1))
#main effects
anova(pgls_photo_bin)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.35757 0.35757 231.3066 3.318e-09 ***
## photo2 2 0.02991 0.01495 9.6738 0.003147 **
## Residuals 12 0.01855 0.00155
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_photo_bin)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + photo2, data = species.comp2,
## lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.084884 -0.027957 -0.003879 0.013107 0.042196
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.10307
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.237405 0.105786 -11.6973 6.425e-08 ***
## log10(length_av) 0.774819 0.071712 10.8045 1.543e-07 ***
## photo2nonpesta 0.104921 0.023898 4.3904 0.0008801 ***
## photo2none 0.050714 0.043692 1.1607 0.2683204
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03932 on 12 degrees of freedom
## Multiple R-squared: 0.9543, Adjusted R-squared: 0.9429
## F-statistic: 83.55 on 3 and 12 DF, p-value: 2.614e-08When lambda is fitted by ML, the estimate for lambda is 0, though this estimate is not significant. So, we fit at both extremes of lambda to see the range of possible results.
Test where lambda = 0
Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_photo0_bin <- pgls(log10(eye_av) ~ log10(length_av) + photo2,
data = species.comp2,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo0_bin)
par(mfrow = c(1, 1))
#main effects
anova(pgls_photo0_bin)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.35757 0.35757 231.3046 3.318e-09 ***
## photo2 2 0.02991 0.01495 9.6737 0.003147 **
## Residuals 12 0.01855 0.00155
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_photo0_bin)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + photo2, data = species.comp2,
## lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.084884 -0.027957 -0.003879 0.013107 0.042196
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.237405 0.105786 -11.6973 6.425e-08 ***
## log10(length_av) 0.774819 0.071712 10.8045 1.543e-07 ***
## photo2nonpesta 0.104921 0.023898 4.3904 0.0008801 ***
## photo2none 0.050714 0.043692 1.1607 0.2683199
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03932 on 12 degrees of freedom
## Multiple R-squared: 0.9543, Adjusted R-squared: 0.9429
## F-statistic: 83.55 on 3 and 12 DF, p-value: 2.614e-08Photophore type has a significant effect on relative eye size when lambda = 0.
Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 1)
pgls_photo1_bin <- pgls(log10(eye_av) ~ log10(length_av) + photo2,
data = species.comp2,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo1_bin)
par(mfrow = c(1, 1))
#main effects
anova(pgls_photo1_bin)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.253779 0.253779 115.4961 1.636e-07 ***
## photo2 2 0.021624 0.010812 4.9206 0.02751 *
## Residuals 12 0.026368 0.002197
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#coefficients and contrasts
summary(pgls_photo1_bin)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + photo2, data = species.comp2,
## lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06902 -0.02289 -0.00536 0.03717 0.06538
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.191766 0.122725 -9.7109 4.912e-07 ***
## log10(length_av) 0.744165 0.080955 9.1923 8.827e-07 ***
## photo2nonpesta 0.106186 0.034696 3.0605 0.00989 **
## photo2none 0.066568 0.047496 1.4016 0.18638
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04688 on 12 degrees of freedom
## Multiple R-squared: 0.9126, Adjusted R-squared: 0.8908
## F-statistic: 41.78 on 3 and 12 DF, p-value: 1.255e-06Here when we bin into 3 groups for photophores (organ of pesta, no organ of pesta, no photophores at all), we find that organ type has a significant effect on eye size in both models (lambda = 1 and lambda = 0). This is a robust result, as whatever the true value of lambda is, the association should hold. The only caveat I will add to this is that since we only have a single evolutionary transition from organ of pesta to non organ of pesta represented in our dataset, the association would be more convincing if found later across larger samples/trees/etc (just something small to mention in discussion).
Depth distributions
Does eye size correlate with daytime or nighttime depth distributions? Are these relationships affected by vertical migration?
Prediction: Eye size increases with daytime depth, as these animals may require the greatest sensitivity to light to cope with the exponential attenuation of light with depth.
Note: Did we have predictions about how vertical migration (which we are factoring in as the interaction effect of day and night depths) might affect eye size? I might predict that stronger migrators might have larger eyes, as they may be able to use ambient light for vision more of the time (things that stay deep at night would not have ambient light available at night). Or they might be using their eyes to follow an isolume. Just a couple ideas!
Exploratory plots (depth)
Before incorporating phylogeny, we can inspect the distribution of our data with scatter plots for absolute and relative eye size and body size vs. daytime and nighttime depth. In the below interactive plots, daytime depths are orange and nighttime depths are blue. Maximum reported depths are in black squares. Depth is not logged.
Absolute eye size with depth.
# Absolute eye size with depths ------
plot_abs_depths <- ggplot(species.phy, aes(x = eye_av, y = Day_Avg, text = genus_species)) +
geom_point(size = 2, alpha = 0.8, col = "orange") +
scale_y_reverse(name = "Depth (m)") +
xlab("Eye diameter (mm)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_point(aes(x = eye_av, y = Night_Avg), size = 2, alpha = 0.8, col = "royalblue2") +
geom_point(aes(x = eye_av, y = Max_reported), size = 1.5, alpha = 0.8, pch = 15, col = "black")
ggplotly(plot_abs_depths)Body size with depth.
# Body size with depth ------
plot_length_depths <- ggplot(species.phy, aes(x = length_av, y = Day_Avg, text = genus_species)) +
geom_point(size = 2, alpha = 0.8, col = "orange") +
scale_y_reverse(name = "Depth (m)") +
xlab("Body length (mm)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_point(aes(x = length_av, y = Night_Avg), size = 2, alpha = 0.8, col = "royalblue2")+
geom_point(aes(x = length_av, y = Max_reported), size = 1.5, alpha = 0.8, pch = 15, col = "black")
ggplotly(plot_length_depths)Relative eye size (eye diameter/body length) with depth.
# Eye size/body size (relative eye size) with max depth -----
plot_rel_depth <- ggplot(species.phy, aes(x = eye_av/length_av, y = Day_Avg, text = genus_species)) +
geom_point(size = 2, alpha = 0.8, col = "orange") +
scale_y_reverse(name = "Depth (m)") +
xlab("Eye/body ratio") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
geom_point(aes(x = eye_av/length_av, y = Night_Avg), size = 2, alpha = 0.8, col = "royalblue2") +
geom_point(aes(x = eye_av/length_av, y = Max_reported), size = 1.5, alpha = 0.8, pch = 15, col = "black")
ggplotly(plot_rel_depth)Is daytime depth correlated with maximum depth? (Presumably these are two measures of a similar thing)
#is daytime depth correlated with maximum depth?
plot_daymax <- ggplot(species.phy, aes(x = Day_Avg, y = Max_reported, text = genus_species)) +
geom_point(size = 2, alpha = 0.8, col = "black") +
scale_y_reverse(name ="Max reported depth (m)") +
xlab("Daytime mean depth (m)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplotly(plot_daymax)Visual inspection of the data shows that absolute eye size, body size, and relative eye size maybe could seem to increase with daytime depth, but if so it’s only driven by a couple of data points. Nighttime depth doesn’t vary much across species, but to me it looks like the species with larger body sizes might also have the larger differences between nighttime and daytime depths. Maximum reported depth shows a similar pattern to daytime depth. We next test this in an evolutionary framework.
Mean depths
Absolute eye size vs. mean day and night depths
Here, we run models of eye size vs. day x night mean depths, with lambda bounded at 0 and then again at 1.
Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_depth0 <- pgls(eye_av ~ Day_Avg * Night_Avg,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depth0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_depth0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Day_Avg 1 0.60507 0.60507 10.6842 0.006720 **
## Night_Avg 1 0.00627 0.00627 0.1108 0.744988
## Day_Avg:Night_Avg 1 0.87460 0.87460 15.4435 0.001999 **
## Residuals 12 0.67958 0.05663
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#coefficients and contrasts
summary(pgls_depth0)##
## Call:
## pgls(formula = eye_av ~ Day_Avg * Night_Avg, data = species.comp,
## lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25976 -0.04300 0.05435 0.19629 0.41490
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.8839e-01 3.0478e-01 -1.2744 0.226663
## Day_Avg 1.0386e-03 2.4323e-04 4.2699 0.001088 **
## Night_Avg 6.2419e-03 1.6506e-03 3.7815 0.002617 **
## Day_Avg:Night_Avg -3.6627e-06 9.3203e-07 -3.9298 0.001999 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.238 on 12 degrees of freedom
## Multiple R-squared: 0.6862, Adjusted R-squared: 0.6077
## F-statistic: 8.746 on 3 and 12 DF, p-value: 0.002394There is a significant effect of daytime depth and a significant interaction effect of day X night depth on absolute eye size when lambda = 0.
Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 1)
pgls_depth1 <- pgls(eye_av ~ Day_Avg * Night_Avg,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depth1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_depth1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Day_Avg 1 0.05746 0.05746 1.0034 0.336259
## Night_Avg 1 0.00178 0.00178 0.0312 0.862811
## Day_Avg:Night_Avg 1 0.97260 0.97260 16.9856 0.001417 **
## Residuals 12 0.68712 0.05726
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#coefficients and contrasts
summary(pgls_depth1)##
## Call:
## pgls(formula = eye_av ~ Day_Avg * Night_Avg, data = species.comp,
## lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.47134 -0.08056 0.01827 0.11141 0.30318
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.0919e-01 2.7596e-01 -0.3957 0.699282
## Day_Avg 6.6304e-04 2.0266e-04 3.2717 0.006682 **
## Night_Avg 5.8104e-03 1.4414e-03 4.0310 0.001666 **
## Day_Avg:Night_Avg -3.3074e-06 8.0249e-07 -4.1214 0.001417 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2393 on 12 degrees of freedom
## Multiple R-squared: 0.6003, Adjusted R-squared: 0.5003
## F-statistic: 6.007 on 3 and 12 DF, p-value: 0.009693When lambda = 1, only the interaction between day and night depths have a significant effect on absolute eye size.
Relative eye size vs. day and night depths.
Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with covariates of daytime depth, nighttime depth, and the interaction between daytime and nighttime depths.
#run PGLS model using the Maximum Liklihood estimate of lambda
pgls_depth <- pgls(log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg,
data = species.comp,
lambda = "ML", bounds = list(lambda = c(1e-05, 1)),
#uses Maximum Liklihood estimate of lambda
param.CI = 0.95)
#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_depth) #plot the linear model
par(mfrow = c(1,1)) #set plotting window back to one plotThe model checks here look a bit weird (residuals look bimodal; assumption is normality).
Next we can take a look at the overall significance of the main effects/interactions and sequential sums of squares. Here we can see that body length has a significant effect on eye size when depth is corrected for (expected, and we already knew that), but daytime depth, nighttime depth, and the interaction between daytime and nighttime depth do not have significant effects on relative eye size.
#anova
anova(pgls_depth)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.253779 0.253779 74.6881 3.112e-06 ***
## Day_Avg 1 0.004141 0.004141 1.2187 0.2932
## Night_Avg 1 0.006218 0.006218 1.8299 0.2033
## Day_Avg:Night_Avg 1 0.000256 0.000256 0.0755 0.7886
## Residuals 11 0.037376 0.003398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Finally, we look at the PGLS coefficients as well significance with respect to contrasts.
#summary
summary(pgls_depth)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg,
## data = species.comp, lambda = "ML", param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.075678 -0.052157 0.001252 0.039596 0.061882
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.19099
## upper bound : 1.000, p = 1
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3424e+00 1.7428e-01 -7.7025 9.348e-06 ***
## log10(length_av) 8.6045e-01 1.5864e-01 5.4240 0.0002089 ***
## Day_Avg 5.1702e-05 6.4331e-05 0.8037 0.4386036
## Night_Avg -2.4008e-04 5.7990e-04 -0.4140 0.6868326
## Day_Avg:Night_Avg 8.7933e-08 3.2009e-07 0.2747 0.7886296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05829 on 11 degrees of freedom
## Multiple R-squared: 0.8761, Adjusted R-squared: 0.8311
## F-statistic: 19.45 on 4 and 11 DF, p-value: 5.969e-05The output also gives a maximum-likelihood estimation of lambda = 1.00, indicating that there is high phylogenetic signal in the residuals of the model (and that caper used a lambda of 1 to set the variance-covariance matrix in the regression). However, small sample sizes have no power to estimate the true value of lambda, so this is not very meaningful in this dataset. We should decide on a strategy for analysis (as mentioned above, may want to run all models again with lambda fixed at 0 as a comparison for the supp info so we can discuss how that affects our results).
#Likelihood plot for Pagel's lambda
plot(pgls.profile(pgls_depth, "lambda"),
main = "Pagel's lambda",
sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")
Bounded extremes of mean depth models
With sample sizes too small to generate a robust estimate of lambda, we here again run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible parameter estimates are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.
Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_depth0 <- pgls(log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depth0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_depth0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.35757 0.35757 113.8476 3.853e-07 ***
## Day_Avg 1 0.00801 0.00801 2.5510 0.1385
## Night_Avg 1 0.00581 0.00581 1.8509 0.2009
## Day_Avg:Night_Avg 1 0.00009 0.00009 0.0272 0.8719
## Residuals 11 0.03455 0.00314
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#coefficients and contrasts
summary(pgls_depth0)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg,
## data = species.comp, lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.086086 -0.012489 0.003432 0.033701 0.096951
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3733e+00 1.5202e-01 -9.0334 2.022e-06 ***
## log10(length_av) 8.4631e-01 1.4585e-01 5.8028 0.0001188 ***
## Day_Avg 1.1361e-04 8.0869e-05 1.4049 0.1876539
## Night_Avg -1.3455e-05 5.7325e-04 -0.0235 0.9816948
## Day_Avg:Night_Avg -5.3441e-08 3.2387e-07 -0.1650 0.8719291
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05604 on 11 degrees of freedom
## Multiple R-squared: 0.9149, Adjusted R-squared: 0.884
## F-statistic: 29.57 on 4 and 11 DF, p-value: 7.848e-06Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals). *Note that when we used a ML estimate of lambda, this was the value used, though the lambda estimate was not significant.
#fit model for eye diameter ~ body length (lambda = 1)
pgls_depth1 <- pgls(log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depth1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_depth1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.253779 0.253779 74.6881 3.112e-06 ***
## Day_Avg 1 0.004141 0.004141 1.2187 0.2932
## Night_Avg 1 0.006218 0.006218 1.8299 0.2033
## Day_Avg:Night_Avg 1 0.000256 0.000256 0.0755 0.7886
## Residuals 11 0.037376 0.003398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_depth1)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg,
## data = species.comp, lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.075678 -0.052157 0.001252 0.039596 0.061882
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3424e+00 1.7428e-01 -7.7025 9.348e-06 ***
## log10(length_av) 8.6045e-01 1.5864e-01 5.4240 0.0002089 ***
## Day_Avg 5.1702e-05 6.4331e-05 0.8037 0.4386036
## Night_Avg -2.4008e-04 5.7990e-04 -0.4140 0.6868326
## Day_Avg:Night_Avg 8.7933e-08 3.2009e-07 0.2747 0.7886296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05829 on 11 degrees of freedom
## Multiple R-squared: 0.8761, Adjusted R-squared: 0.8311
## F-statistic: 19.45 on 4 and 11 DF, p-value: 5.969e-05In both models bounding the possible true value of lambda, daytime depth, nighttime depth, and the interaction between day and night depths do not have significant effects on relative eye size in these species. This means that our finding is robust, regardless of our challenge estimating phylogenetic signal in our model residuals.
Maximum depth
While we found no association between eye size and mean daytime or nighttime depths, we can also test eye size vs. the maximum reported depth for each species.
Absolute eye size vs. max depth
Here is a model of eye size ~ maximum reported depth fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_maxdepth0 <- pgls(eye_av ~ Max_reported,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_maxdepth0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_maxdepth0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Max_reported 1 1.0314 1.03141 12.732 0.003087 **
## Residuals 14 1.1341 0.08101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#coefficients and contrasts
summary(pgls_maxdepth0)##
## Call:
## pgls(formula = eye_av ~ Max_reported, data = species.comp, lambda = 1e-05,
## param.CI = 0.95, bounds = list(lambda = c(1e-05, 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35034 -0.20774 -0.03950 0.05696 0.58971
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.26375005 0.23390646 1.1276 0.278449
## Max_reported 0.00049307 0.00013818 3.5682 0.003087 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2846 on 14 degrees of freedom
## Multiple R-squared: 0.4763, Adjusted R-squared: 0.4389
## F-statistic: 12.73 on 1 and 14 DF, p-value: 0.003087Maximum reported depth has a significant effect on absolute eye size when lambda = 0.
Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 1)
pgls_maxdepth1 <- pgls(eye_av ~ Max_reported,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_maxdepth1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_maxdepth1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Max_reported 1 0.45511 0.45511 5.0414 0.04142 *
## Residuals 14 1.26385 0.09028
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_maxdepth1)##
## Call:
## pgls(formula = eye_av ~ Max_reported, data = species.comp, lambda = 1,
## param.CI = 0.95, bounds = list(lambda = c(1, 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.66883 -0.28878 -0.07065 0.02554 0.37075
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.42595485 0.27755652 1.5347 0.14715
## Max_reported 0.00036789 0.00016385 2.2453 0.04142 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3005 on 14 degrees of freedom
## Multiple R-squared: 0.2648, Adjusted R-squared: 0.2122
## F-statistic: 5.041 on 1 and 14 DF, p-value: 0.04142Maximum reported depth has a significant effect on absolute eye size when lambda = 1.
In both models bounding the possible true value of lambda, maximum depth has a significant effect on absolute eye size.
Relative eye size vs. max depth
Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with covariate of maximum depth.
#run PGLS model using the Maximum Liklihood estimate of lambda
pgls_maxdepth <- pgls(log10(eye_av) ~ log10(length_av) + Max_reported,
data = species.comp,
lambda = "ML", bounds = list(lambda = c(1e-05, 1)),
#uses Maximum Liklihood estimate of lambda
param.CI = 0.95)
#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_maxdepth) #plot the linear model
par(mfrow = c(1,1)) #set plotting window back to one plotThe model checks here look ok.
Next we can take a look at the overall significance of the main effects/interactions and sequential sums of squares. Here we can see that body length has a significant effect on eye size when maximum reported depth is corrected for (expected, and we already knew that), but that maximum depth does not have a significant effect on eye size.
#anova
anova(pgls_maxdepth)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.89, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.257318 0.257318 82.669 5.364e-07 ***
## Max_reported 1 0.006247 0.006247 2.007 0.1801
## Residuals 13 0.040464 0.003113
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Finally, we look at the PGLS coefficients as well significance with respect to contrasts.
#summary
summary(pgls_maxdepth)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Max_reported,
## data = species.comp, lambda = "ML", param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.05638 -0.02305 0.01184 0.06568 0.08418
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.886
## lower bound : 0.000, p = 0.28682
## upper bound : 1.000, p = 0.80796
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.2799e+00 1.4323e-01 -8.9360 6.528e-07 ***
## log10(length_av) 7.6957e-01 1.0295e-01 7.4751 4.663e-06 ***
## Max_reported 4.8325e-05 3.4111e-05 1.4167 0.1801
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05579 on 13 degrees of freedom
## Multiple R-squared: 0.8669, Adjusted R-squared: 0.8464
## F-statistic: 42.34 on 2 and 13 DF, p-value: 2.028e-06The output also gives a maximum-likelihood estimation of lambda = 0.89, indicating that there is high phylogenetic signal in the residuals of the model (and that caper used a lambda of 0.89 to set the variance-covariance matrix in the regression). However, small sample sizes have no power to estimate the true value of lambda, so this is not very meaningful in this dataset. We additionally run models with lambda fixed at 0 and at 1 as comparisons for the supp info so we can discuss how that affects our results.
#Likelihood plot for Pagel's lambda
plot(pgls.profile(pgls_maxdepth, "lambda"),
main = "Pagel's lambda",
sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")
Bounded extremes of maximum reported depth models
With sample sizes too small to generate a robust estimate of lambda, we here run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible results are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.
Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_maxdepth0 <- pgls(log10(eye_av) ~ log10(length_av) + Max_reported,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_maxdepth0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_maxdepth0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.35757 0.35757 123.7525 5.157e-08 ***
## Max_reported 1 0.01090 0.01090 3.7715 0.07412 .
## Residuals 13 0.03756 0.00289
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_maxdepth0)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Max_reported,
## data = species.comp, lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.09807 -0.03771 0.01952 0.04081 0.07693
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3396e+00 1.3743e-01 -9.7469 2.421e-07 ***
## log10(length_av) 8.0134e-01 1.0607e-01 7.5547 4.162e-06 ***
## Max_reported 6.4411e-05 3.3167e-05 1.9420 0.07412 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05375 on 13 degrees of freedom
## Multiple R-squared: 0.9075, Adjusted R-squared: 0.8933
## F-statistic: 63.76 on 2 and 13 DF, p-value: 1.907e-07When lambda = 0, maximum depth has no effect on relative eye size.
Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 1)
pgls_maxdepth1 <- pgls(log10(eye_av) ~ log10(length_av) + Max_reported,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_maxdepth1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_maxdepth1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.253779 0.253779 78.0290 7.436e-07 ***
## Max_reported 1 0.005711 0.005711 1.7559 0.208
## Residuals 13 0.042281 0.003252
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_maxdepth1)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Max_reported,
## data = species.comp, lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.072284 -0.035744 0.008099 0.045397 0.078221
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.2672e+00 1.4497e-01 -8.7411 8.369e-07 ***
## log10(length_av) 7.6322e-01 1.0345e-01 7.3773 5.367e-06 ***
## Max_reported 4.5769e-05 3.4539e-05 1.3251 0.208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05703 on 13 degrees of freedom
## Multiple R-squared: 0.8599, Adjusted R-squared: 0.8383
## F-statistic: 39.89 on 2 and 13 DF, p-value: 2.832e-06When lambda = 1, maximum depth has no effect on relative eye size.
In both models bounding the possible true value of lambda, maximum depth does not have significant effects on relative eye size in these species.
Depth range traversed in DVM (daytime mean - nighttime mean)
While we found no association between eye size and mean daytime or nighttime depths or the interaction between them, I also wanted to check whether there was an association between eye size and the distance traversed during DVM.
Absolute eye size vs. depths traversed
Here is the model of eye size vs. depth range fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_depthrange0 <- pgls(eye_av ~ Range,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_depthrange0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Range 1 0.4582 0.45820 3.7572 0.07301 .
## Residuals 14 1.7073 0.12195
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_depthrange0)##
## Call:
## pgls(formula = eye_av ~ Range, data = species.comp, lambda = 1e-05,
## param.CI = 0.95, bounds = list(lambda = c(1e-05, 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62839 -0.26817 0.04641 0.22210 0.56496
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.78225824 0.16726887 4.6767 0.0003566 ***
## Range 0.00055569 0.00028668 1.9384 0.0730117 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3492 on 14 degrees of freedom
## Multiple R-squared: 0.2116, Adjusted R-squared: 0.1553
## F-statistic: 3.757 on 1 and 14 DF, p-value: 0.07301When lambda = 0, depth range has no significant effect on absolute eye size.
Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 1)
pgls_depthrange1 <- pgls(eye_av ~ Range,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_depthrange1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Range 1 0.04454 0.044539 0.3724 0.5515
## Residuals 14 1.67442 0.119602#parameter estimates
summary(pgls_depthrange1)##
## Call:
## pgls(formula = eye_av ~ Range, data = species.comp, lambda = 1,
## param.CI = 0.95, bounds = list(lambda = c(1, 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64588 -0.44759 -0.13224 0.03535 0.27527
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.93155235 0.16629239 5.6019 6.524e-05 ***
## Range 0.00013888 0.00022758 0.6102 0.5515
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3458 on 14 degrees of freedom
## Multiple R-squared: 0.02591, Adjusted R-squared: -0.04367
## F-statistic: 0.3724 on 1 and 14 DF, p-value: 0.5515When lambda = 0, depth range has no significant effect on absolute eye size. Both models agree there is no effect of depth traversed during migration on absolute eye size.
Relative eye size vs. depth range traversed
Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with depth range traversed during DVM as a covariate.
#run PGLS model using the Maximum Liklihood estimate of lambda
pgls_depthrange <- pgls(log10(eye_av) ~ log10(length_av) + Range,
data = species.comp,
lambda = "ML", bounds = list(lambda = c(1e-05, 1)),
#uses Maximum Liklihood estimate of lambda
param.CI = 0.95)
#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_depthrange) #plot the linear model
par(mfrow = c(1,1)) #set plotting window back to one plotThe model checks here look a bit wonky.
Next we can take a look at the overall significance of the main effects and sequential sums of squares. Here we can see that body length has a significant effect on eye size when depth range is corrected for (expected, and we already knew that), but that depth range does not have a significant effect on eye size.
#anova
anova(pgls_depthrange)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.253779 0.253779 86.7481 4.078e-07 ***
## Range 1 0.009961 0.009961 3.4048 0.0879 .
## Residuals 13 0.038031 0.002925
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Finally, we look at the PGLS coefficients as well significance with respect to contrasts.
#summary
summary(pgls_depthrange)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Range, data = species.comp,
## lambda = "ML", param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.069275 -0.051016 0.008185 0.041543 0.058192
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.22392
## upper bound : 1.000, p = 1
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3117e+00 1.3728e-01 -9.5547 3.044e-07 ***
## log10(length_av) 8.1984e-01 8.8362e-02 9.2782 4.261e-07 ***
## Range 6.5688e-05 3.5599e-05 1.8452 0.0879 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05409 on 13 degrees of freedom
## Multiple R-squared: 0.874, Adjusted R-squared: 0.8546
## F-statistic: 45.08 on 2 and 13 DF, p-value: 1.422e-06The output also gives a maximum-likelihood estimation of lambda = 1, indicating that there is high phylogenetic signal in the residuals of the model (and that caper used a lambda of 1 to set the variance-covariance matrix in the regression). However, small sample sizes have no power to estimate the true value of lambda, so this is not very meaningful in this dataset. We additionally run models with lambda fixed at 0 and at 1 as comparisons for the supp info so we can discuss how that affects our results.
#Likelihood plot for Pagel's lambda
plot(pgls.profile(pgls_depthrange, "lambda"),
main = "Pagel's lambda",
sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")
Bounded extremes of depth range models
With sample sizes too small to generate a robust estimate of lambda, we here run models with the lowest and highest possible bounds of lambda (0 and 1), so that we can see what the range of possible results are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.
Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_depthrange0 <- pgls(log10(eye_av) ~ log10(length_av) + Range,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_depthrange0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.35757 0.35757 134.153 3.192e-08 ***
## Range 1 0.01381 0.01381 5.181 0.0404 *
## Residuals 13 0.03465 0.00267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_depthrange0)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Range, data = species.comp,
## lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.088931 -0.012322 0.002289 0.031130 0.099003
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3794e+00 1.2604e-01 -10.9445 6.253e-08 ***
## log10(length_av) 8.6156e-01 8.5383e-02 10.0905 1.621e-07 ***
## Range 1.0275e-04 4.5143e-05 2.2762 0.0404 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05163 on 13 degrees of freedom
## Multiple R-squared: 0.9147, Adjusted R-squared: 0.9015
## F-statistic: 69.67 on 2 and 13 DF, p-value: 1.128e-07When lambda = 0, depth range traversed has a significant effect on relative eye size.
Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 1)
pgls_depthrange1 <- pgls(log10(eye_av) ~ log10(length_av) + Range,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_depthrange1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.253779 0.253779 86.7481 4.078e-07 ***
## Range 1 0.009961 0.009961 3.4048 0.0879 .
## Residuals 13 0.038031 0.002925
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_depthrange1)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Range, data = species.comp,
## lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.069275 -0.051016 0.008185 0.041543 0.058192
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3117e+00 1.3728e-01 -9.5547 3.044e-07 ***
## log10(length_av) 8.1984e-01 8.8362e-02 9.2782 4.261e-07 ***
## Range 6.5688e-05 3.5599e-05 1.8452 0.0879 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05409 on 13 degrees of freedom
## Multiple R-squared: 0.874, Adjusted R-squared: 0.8546
## F-statistic: 45.08 on 2 and 13 DF, p-value: 1.422e-06In the model with lambda set to 0, depth range does not have a significant effect on relative eye size. Thus, the answer depends on the true value of lambda (which we can’t get at with our small sample size).
Let’s just see what happens with a lambda of 0.5 for better understanding.
#fit model for eye diameter ~ body length (lambda = 1)
pgls_depthrange05 <- pgls(log10(eye_av) ~ log10(length_av) + Range,
data = species.comp,
lambda = 0.5, #set lambda to 0.5
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_depthrange05)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.50, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.280831 0.280831 103.4663 1.476e-07 ***
## Range 1 0.010073 0.010073 3.7112 0.0762 .
## Residuals 13 0.035285 0.002714
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_depthrange05)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Range, data = species.comp,
## lambda = 0.5, param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06357 -0.02033 0.03424 0.05163 0.08144
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.500
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3526e+00 1.3250e-01 -10.2085 1.415e-07 ***
## log10(length_av) 8.4502e-01 8.7069e-02 9.7052 2.543e-07 ***
## Range 8.2101e-05 4.2618e-05 1.9264 0.0762 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0521 on 13 degrees of freedom
## Multiple R-squared: 0.8918, Adjusted R-squared: 0.8752
## F-statistic: 53.59 on 2 and 13 DF, p-value: 5.27e-07At lambda = 0.5 depth range does not have a significant effect on relative eye size.
Habitat
Preliminary plots (habitat)
Before incorporating phylogeny, we can inspect the distribution of our data with boxplots for absolute and relative eye size be habitat.
# Phylogeny with eye size and investment and ecology -------
#color vector for photophores
col_eco <- c("no" = "#7FDBFF",
"yes" = "#FF851B")
#resort trait dataset to the order of tree tip labels
species.phy <- species.phy[tree.plot$tip.label, ]
#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye,
scale = 1.5,
tip.labels = TRUE,
col = unname(col_eco[as.factor(species.phy$Benthopelagic)]),
offset = 1,
ftype = "bi",
fsize = 0.7)
#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
pch = 19, #shape of labels
col = "mediumpurple",
offset = 0)
#add legend for photophores
legend(x = "bottomleft", legend = c("pelagic", "benthopelagic"), pch = 22, pt.cex= 2, pt.bg = col_eco, cex = 1, bty = "n", horiz = F)
#Boxplots for absolute eye size across photophores -----
plot_abs_eco <- ggplot(data = species.phy, aes(x = reorder(Benthopelagic, eye_av, FUN = mean), y = eye_av)) +
geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
geom_jitter(shape = 19, size = 1.3, alpha = 0.5, position = position_jitter(0.15)) +
theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black")) + #controls background
xlab("Benthopelagic") +
ylab("Eye diameter (mm)") +
ggtitle("Absolute eye size")
ggplotly(plot_abs_eco)#Boxplots for relative eye investment across photophores -----
plot_rel_eco <- ggplot(data = species.phy, aes(x = reorder(Benthopelagic, pglsres, FUN = mean), y = pglsres)) +
geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
geom_jitter(shape = 19, size = 1.3, alpha = 0.5, position = position_jitter(0.15)) +
theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black")) + #controls background
xlab("Benthopelagic") +
ylab("Eye ~ body residuals") +
ggtitle("Relative eye investment")
ggplotly(plot_rel_eco) #Export figure------
pdf("../Figures/phylo_eco.pdf", width = 6, height = 4)
#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye,
scale = 1.5,
tip.labels = TRUE,
col = unname(col_eco[as.factor(species.phy$Benthopelagic)]),
offset = 1,
ftype = "bi",
fsize = 0.7)
#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
pch = 19, #shape of labels
col = "mediumpurple",
offset = 0)
#add legend for ecology
legend(x = "bottomleft", legend = c("pelagic", "benthopelagic"), pch = 22, pt.cex= 2, pt.bg = col_eco, cex = 1, bty = "n", horiz = F)
dev.off()It looks like benthopelagic species have big-ish eyes, but it’s small sample size so hard to say. Distribution overlaps with pelagic species. We next test this in an evolutionary framework.
Absolute eye size vs. habitat
We test for effects of habitat on absolute eye size using a PGLS of eye size ~ habitat.
Here is a model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_eco0 <- pgls(eye_av ~ as.factor(Benthopelagic),
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_eco0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_eco0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Benthopelagic) 1 0.11959 0.11959 0.8183 0.381
## Residuals 14 2.04594 0.14614#parameter estimates
summary(pgls_eco0)##
## Call:
## pgls(formula = eye_av ~ as.factor(Benthopelagic), data = species.comp,
## lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59327 -0.29244 0.02655 0.26145 0.58982
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.01729 0.10603 9.5947 1.556e-07 ***
## as.factor(Benthopelagic)yes 0.22150 0.24486 0.9046 0.381
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3823 on 14 degrees of freedom
## Multiple R-squared: 0.05522, Adjusted R-squared: -0.01226
## F-statistic: 0.8183 on 1 and 14 DF, p-value: 0.381When lambda = 0, habitat has no effect on absolute eye size.
Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 1)
pgls_eco1 <- pgls(eye_av ~ as.factor(Benthopelagic),
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_eco1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_eco1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: eye_av
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Benthopelagic) 1 0.00128 0.001281 0.0104 0.9201
## Residuals 14 1.71768 0.122692#parameter estimates
summary(pgls_eco1)##
## Call:
## pgls(formula = eye_av ~ as.factor(Benthopelagic), data = species.comp,
## lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71189 -0.44229 -0.13957 0.03058 0.32764
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.988406 0.138957 7.1130 5.229e-06 ***
## as.factor(Benthopelagic)yes 0.021808 0.213436 0.1022 0.9201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3503 on 14 degrees of freedom
## Multiple R-squared: 0.0007451, Adjusted R-squared: -0.07063
## F-statistic: 0.01044 on 1 and 14 DF, p-value: 0.9201Habitat (benthopelagic/pelagic) does not have a significant effect on absolute eye size when lambda is set to 0 or 1. Thus, our finding is robust whatever the true value of lambda is. However, with n = 3 for benthopelagic species, we should note that our power here is very low.
Relative eye size vs. habitat
Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with habitat (benthopelagic vs. pelagic) as a covariate.
#run PGLS modelusing the Maximum Liklihood estimate of lambda
pgls_eco <- pgls(log10(eye_av) ~ log10(length_av) + Benthopelagic,
data = species.comp,
lambda = "ML",
#uses Maximum Liklihood estimate of lambda
param.CI = 0.95)
#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_eco) #plot the linear model
par(mfrow = c(1,1)) #set plotting window back to one plotNext we can take a look at the overall significance of the main, 2-way, 3-way, etc. interactions and sequential sums of squares. Here we can see that body length has a significant effect on eye size (expected, we already knew that), but ecology (benthopelagic or not) has no effect on eye size (p = 0.57).
#anova
anova(pgls_eco)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.253779 0.253779 70.4597 1.316e-06 ***
## Benthopelagic 1 0.001169 0.001169 0.3245 0.5786
## Residuals 13 0.046823 0.003602
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Finally, we look at the PGLS coefficients as well as significance with respect to contrasts.
#summary
summary(pgls_eco)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Benthopelagic,
## data = species.comp, lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.092420 -0.049981 0.009349 0.043498 0.080494
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.11794
## upper bound : 1.000, p = 1
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.287737 0.151661 -8.4909 1.158e-06 ***
## log10(length_av) 0.820905 0.098088 8.3691 1.360e-06 ***
## Benthopelagicyes 0.020845 0.036592 0.5697 0.5786
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06001 on 13 degrees of freedom
## Multiple R-squared: 0.8448, Adjusted R-squared: 0.821
## F-statistic: 35.39 on 2 and 13 DF, p-value: 5.496e-06The output also includes the maximum-likelihood estimation of lambda, which represents the phylogenetic signal for PGLS model residuals. Here, the estimate for lambda is 1.00, which means that there is the highest phylogenetic signal. However, as mentioned above, when you have <20 species the power to estimate lambda is not there, so below we will run models with the min and max possible values of lambda to see how that affects our results.
#Likelihood plot for Pagel's lambda
plot(pgls.profile(pgls_eco, "lambda"),
main = "Pagel's lambda",
sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")
Bounded extremes of habitat models
With sample sizes too small to generate a robust estimate of lambda, one option is to run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible parameter estimates are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.
Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).
#fit model for eye diameter ~ body length (lambda = 0)
pgls_eco0 <- pgls(log10(eye_av) ~ log10(length_av) + Benthopelagic,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_eco0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_eco0)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.35757 0.35757 102.5837 1.551e-07 ***
## Benthopelagic 1 0.00315 0.00315 0.9026 0.3594
## Residuals 13 0.04531 0.00349
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_eco0)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Benthopelagic,
## data = species.comp, lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05,
## 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.102084 -0.026881 0.005142 0.048346 0.095391
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.409336 0.143626 -9.8125 2.24e-07 ***
## log10(length_av) 0.909764 0.093762 9.7029 2.55e-07 ***
## Benthopelagicyes 0.036746 0.038678 0.9501 0.3594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05904 on 13 degrees of freedom
## Multiple R-squared: 0.8884, Adjusted R-squared: 0.8712
## F-statistic: 51.74 on 2 and 13 DF, p-value: 6.454e-07With lambda = 0, habitat has no effect on relative eye size.
Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals). *Note that when we used a ML estimate of lambda, this was the value used, though the lambda estimate was not significant.
#fit model for eye diameter ~ body length (lambda = 1)
pgls_eco1 <- pgls(log10(eye_av) ~ log10(length_av) + Benthopelagic,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_eco1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_eco1)## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: log10(eye_av)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(length_av) 1 0.253779 0.253779 70.4597 1.316e-06 ***
## Benthopelagic 1 0.001169 0.001169 0.3245 0.5786
## Residuals 13 0.046823 0.003602
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#parameter estimates
summary(pgls_eco1)##
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Benthopelagic,
## data = species.comp, lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1,
## 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.092420 -0.049981 0.009349 0.043498 0.080494
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.287737 0.151661 -8.4909 1.158e-06 ***
## log10(length_av) 0.820905 0.098088 8.3691 1.360e-06 ***
## Benthopelagicyes 0.020845 0.036592 0.5697 0.5786
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06001 on 13 degrees of freedom
## Multiple R-squared: 0.8448, Adjusted R-squared: 0.821
## F-statistic: 35.39 on 2 and 13 DF, p-value: 5.496e-06When lambda = 1, habitat has no effect on relative eye size.
Habitat (benthopelagic/pelagic) does not have a significant effect on relative eye size when lambda is set to 0 or 1. Thus, our finding is robust whatever the true value of lambda is. However, note that we only had n = 3 benthopelagic species in our dataset.